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Bh INTRODUCTION 

The procedure for the design of lowpass, highpass. 
bandpass. and bandstop a2 Cage i filters from lowpass 
Misocetyoes if well esteblished, and is based on transforming 
ere Bact oO pigee fr eee pias: response ES meet desired 
Siac Cf ETI Setics. eee oor {aon Can he accomplished in 
either the analoaq (s) domain or the digital (7) domain 


(Fagure 1). 














| 

| llowpass _ | analog | | 
] feel od L—— (ANALOG TRANSFORMATIONS) =) filter | | 
i [prototype ! rr | 
| (J (s domain) | 
| | | | 
| nee ee | 
L — — BILINEAR Le-----—--—-+—--—-— |] RILINEAP oS 
) I'TRANSFORMATION | | TRANSFORMATION | | 
| 
| (z domain) | 
| ay -—*-—- | 
| | homeoss | final ! | 
| magi tal bo (DIGITAL TRANSFORMATIONS) Spl Saeed 1. | 
| | prototype! | filter | | 
| ee | 

Figure 1. Two methods for digital filter desian. 


Note that each of the horizontal arrows actually 
weoeeeest FONT C2 feeremt transtormations:; 

i* Pewsass woroetotyae ts lowness filter 
lowpeses prototype te highpass filter 
lowpass prototype to bandnass filter 
lowpass prototype to bandstop filter 


iS fy ht) +t 


Therefore. the two classical design techniques. analoac- 
Gicital and dicgcital-didgital design, are based on doing the 
freauencyv response transformations in their respective analoa 
or Gigital domains. Since the theory for the standard filter 
types (Butterworth, Chebyshev, and elliptic) is developed in 
the analog or ‘'s' domain, the two methods for designing 
digitai filters differ basically in the sequence of the 
desicn steps. In the analoag-digital method the analog 
prototype 1s transformed and then digitized, while in the 
Gdigitaei-digital method the anaiog prototype 1s digitized and 
then transformed. 

Each of the transformations (arrows in Figure 1) invogieee 


substituting aptunc tien sot ee ox Poe for the appropriate 


variable in the fiiter transfer function. The purposemes 
this thesis is to deveiop efficient alaqorithms to carryeime 
these transformations and then to implement these algorithms 
in a computer program to aid in the design of di@iieae 
filters. 

Computer-aided desian of ‘SdicitaL filters offers a 
Slonificant advantage since it allows the digital filter 
coefficients to be calculated quickly and accurately. (i 
accuracy (Of Sethe ws ae. coefficients 158 an especially 
Slonificant point since the frequency response of the final 
filter is very sensitive to small inaccuracies 1m 


calculated coefficients. As an example. compare the 


frequency response for the filters in Figure 2 and 3. Fiovure 


RO 


2 shows the frequency response magnitude for a 1/2-dB ripple 


8 


Chebyshev bandpass filter, and Figure 3 shows the frequency 


response when one of the coefficients has been changed by 


1/1000th. 








1.80 | 
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FILTER OUTPUT MAGNITUDE ~ 
0.26 0.60 0.75 1.00 126 


0.00 





0.00 0.05 0.10 0.15 O20 0.25 6.30 038 06.40 045 6.50 
DIGITAL FREQUENCY (F/20KHZ) 


Figure 2. Chebyshev Bandpass Filter 


ae hy) 


ro 
KC 


se OZ 


en eee eee ere Uren 2 outs? +o. 4627 7z*-.6746942+.67529 


Picuner 3. 


0.80 0.76 1.00 1.26 1.60 


FILTER OUTPUT MAGNITUDE 


0.28 


i 
7° 
| 


0.00 





06.060 6.066 6.10 O16 OF0 O76 6.350 036 840 6.45. 8.56 
DIGITAL FREQUENCY (F/20KHZ) «i, ! 


Same Filter With Denominator Coefficient of 
Zz changed to -.67594 
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B. GENERAL, DiSCUSs Lon 
Each of the transformations involves Substztucing a nee 
of two polynomials of s or 2 into the faiter trance 


function which is also a ratio of two polynomials Im “sS 61 








namely 
Yee 
Hly)= == p(x) 
Uly) ly= (1 2s) 
a(x) 
The variables ‘'x' and ‘y' can be either ‘s'’ or ‘2' depend 
on the transformation to be performed. Therefore, we have 
either 
@€o + ary + azv* += ec © 
Hdy) = — 2 eeee 
bo + biy + bey? + ... + bnyt | p(x) (1 5e 
Vm 
q(x) 
on 
1S hose) (p(x) )? {(po(x))® 
@€o + ai + a2 See ee ee 
Gy (3a) (Gate) )™ 
(a(x) ) 
H(x)= WWW SS Se ee P (1 ae 
p(x) (p(x) )?4 (ome) )” 
Dow bn + bi eS +. YS eee 
a(x) (a(x) )? (erat)) 
In general, the order of the numerator ‘'m', is less than the 
order of the denominator 'n' (except for the transformations 


in the digital domain where the result of the bilamee 


transformation is to make n=m). We can qeneralize the 


results, however. by assuming that the order or the numerator 
polynomial Y(x) is equal to the order of the denominator 
polynomiai U(x) and makina the additional coefficients 
Qm+i1,am+2,--.--a@n eGual to zero. Then multiplying both the 
numerator and denominator of the above equation by (a(x))® 
gives 

ages 2 an moix) (atx) )o8-* +m. ww. + an (p(x) ) 


a (1.4) 
Dereeeeteeee ~biep (se(ieuiotieimer ) ate . ete bneol p(x) ) 


Same in a different form. 


S Gus Diese) co (au) 9 =! 
SS SS : ie, 52) 


mes Dia)? g(x)" ! 


Mewewtteige the Humerator and denominator of the same form, 
we can use the same transformation for both the numerator and 


Eenecmimnatcor polynomtals of the filter transfer function. 


namely, 
(1.6) 
1=n 
ee ones ai Xt ciex’ tue . .+ Cie’ Sh a= Saar. plx)'aqt(x)s-! 
1=0 
tm TeGuation 1.6. the a are the coefficients of the old 
numerator or denominator polynomials and ai‘ are the new. 


transformed numerator or denominator polynomial coefficients. 
Since the new coefficients ai'., are a linear combination of 


5 


the old coefficients ai. ama the Seckileiemeenor pita 


ao(x), this suggests the use of a matrix tranetermauionm, 
av = aA) cl Je 
Witernc 
a is the coefficient row vector for the old 


numerator or denominator polynomial of 
lenath=n+1. 


a= [an a: a2... @anj 


a ' is the coefficient row vector for the new 
numerator or denominator polynomial of 
leneth=n +2< 


a’ = [a’o ~ae4 (aoc ee ane 


A(n) is the transformation matrix with dimensions 
(nt+1)x(n'+1). The elements of A(n) depend upon the 
coefficients of p(x) and Gt). "Nete that Tien 
rows and columns of A(n) are numbered 0 to mn tien 


BSE ai A(n)i i (12s 


(where A(n)ii denotes the element in row i, 
column 9 Gf Mate wae 


nh is the order of the transfer function before ie 
transformaticne 


1s the order of the transfer funetion after wee 
transformation. Note that n' equals n times the 
highest order of pix) andwane = 


ey 


For exampie. when n=1 and n'‘'=2, p(x) and q(x) are second 


order and Ecuation Ie) yeecomnec. 


lar a A 2 nn eee anc, A(2)a2 | — at 

— — | | = JaO' al’ a2'| 
IAC 2): 6 ee) Ai2):2 | Les <a 
ae = (15h 


ipaaemeace  Onmtne lowpass Gigital. highpass digital. and 
Patiear tract ermations, the two functions p(x) and a(x) are 
Mmirst’order”’ Professor P.H. Moose of the Naval Postgraduate 
School of Monterey has developed the matrix A(n) for these 
two cases by expanding Equation 1.6 using the binomial 
theorem (Reference 1). This method cannot be used for 
transformations involving higher order polynomials such as 
the bandpass and bandstop transformations and it is difficult 
memcdevelop because intricate manipulation of the binomial 
coefficients are involved. The binomial expansion approach. 
therefore, was not pursued. 

A general method has been found that develops the 
meeurred transformation matrices in an iterative manner. In 
other words, the Matrix used to transform a numerator or 
G@enoninator filter polynomial of order n is developed fron 
Memo trts Lor the polynomial of order n-l. This approach is 
eéesy to implement in a computer program and can be used for 
any or the polynomial substitutions needed. The next chapter 
is concerned with developing the transformation matrices for 
each specific transformation. However, aS an introduction to 
tne method used, inet sx the generic case where both p(x) 
ame tx) are second-order polynomials. Since all of the 
analog and digital transformations involve substitutions with 
polynomiais of second order or less, this second-order 
gcqeneric case will be all that is needed to develop the 


specific transformations later. 


With» p(x)+ax?2 &+ bx + c , Gem = dx coe. faa 
substituting into Equation 1.6 for the [imetuwene@er cas cues 
n=1 (remembering that with p(x) and g(x) of second order, the 


resultina polynomial will be second order), we get 


a'(x) = ao' + ai'xX + a2'x? = an (Ax? + ex + f) + a: (ax? +bx+c) 


(1 ie 


Usina Equation 1.8 or 1.9 it 18 easy to pick out the sau 


coefficients A(1); ; (1=07 1 Mee ©. 292) for the first Ogee 
generic case. The coefficients are, 

— coe 

Pe = d | (1 . date 
Ac b a | 

a. — 


Now with n=2 (n'=4) substituting into Equation 1.6 Gale 


(x) = ao’ + a1 'X + a2'xR* + a3 'xX*® + ag 'x4 (7 


= ao (dx+text+fi)*? + ar (ax*+bxt+c) (ax? textf) t+az2 (ax* +bxt+c) 2 


ana the matrix Aes. 


| £2 2ef 2df+e? 2de a? | 

a2 — (fe ect+bf cd+fateb dbt+ea da! ( Lois) 
ee Ze Zactb* Zab a | 
Le 2 


To see how this matrix is derived from the previous one, 
consider the first row of the matrices in Equations 1.2] 
1. Se For both matrices, this row is multiplied by the 4ao 


coefficient in the old “eeiyvnenaa For n=1 we have aop(xXi 


E 


Saemwe Meme nmr omc ithe coefficients of the powers of x (x!) to 
qet the matrix elements A(n)o:;. For n=2, we have ao (p(x) )+¢ 
and we do the same thing only now we have the additional 
mrecror Of p(x). The matrix elements can be calculated by 
hand in this fashion but it requires a lot of algebra. A 
Simpler method is to use the matrix elements of the previous 
matrix and get the appropriate elements by convolving 
coefficients. For example, consider the sequence of elements 


memenow © (the top mowk of Adl). namely 
e e d ella 


To get the elements for row 0 of A(2). we need to multiply 
(convoive mee ficients) by another p(x). Therefore, with 
rows numbered C to n., and columns number 0 to n', we have the 
men ewinc. 


mor mew, 0 Column O, the appropriate matrix element is 


Fae? )o cmm= A e = fe 


e qd (arm Bs) 


rh kh 


For row 0 column 1. the appropriate matrix element is 


Al2)o, = d e = z2fe 
£ e qd ea 6 


rh 


For row 0 column 2, the appropriate matrix element is 


A(2)o2 = a4 G e if = 2fd+e* 
if e d lee tna -) 


For row 0 column 3. the appropriate matrix een 
A(2)o3 = d & if = Zale 
if 2 a (1.28) 


For row 0 column 4. the appropriate matrix element is 


Al2)oa = d & f = qd? 

if e a (1 [ee 
This convolution of the coefficients of the old matrix teen 
the coefficients of p(x) can be carried out for every row of 
the old matrix to obtain all the rows of the new matrix 
except the last row. In the case of the last row, we are 
multiplyina by an additional a(x) rather than an additional 
p(x) and so must obtain the last row of the new matrix by 
COnVOLVInae tne comeHtes ent of the last row of the old matrix 
with the coefficients o1eciw. 


Im Summary we have, 


(3) for rowel to Trowenae.. 


row i of A(n-1) * coefficients of p(x) 


—-row 1 Of ANG! 


(1 2208 
(11) £or rows: 
-row nm of A(n) = row n-1 of A(n-1) * coefficients Of @eamee 
(1 Zam 


where * indicates linear convolution of the coeftfticireneas 
To see how this method of convolving coefficients works, 


consider multiplvinoe two polynomials: 


fli) ="f2( eee (1 (22 
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Seance 8 15 eme GB@udependent variable, no matter what letter 
we use, the polynomials will remain the same. Substituting 


' 


a = tor x in the above equation gives us, 


Piece — Lavz + )f£3(2-' ) e273") 


eplerewe  —'rlow have the familiar form where we can either 


miler ply the “tworfunctions ~of ‘z'~er convolve the sequences 


found by taking the inverse z-transform. Bor’ sf20z-* ) and 
fexwz !?) we fina 

f2(z-1) = ao + a:2-! + azz-?2 + . . . + anzz72 

tong em = bo + by Z-1 + bez-2 + .. . + bez ® IL les 


Pm icmmmaw clear that we can find the coefficients of the 
product f1(x) by convolving the sequence fan] with the 


sequence {bm}! to produce the sequence {ca! where. 


fan} = !@,@:,@,as3, . .. an] 
Mon; = {bo bi ,b2 , Dam. aes . ba | 

and with q=m+on (Clee. 5) 
Wen = (enc) .Co,¢2, . . VY iq) 


Then bv takino the z-transform of the sequence {ca} and 


SlDsvertutmng x' for ‘z= ' we get. 


PS — ere t cy RK + . wt. UFC hq XT (ieee) 


let 


We now have a simple straiqht-forward method of 
developing the necessary transformation matrices for digital 
fiiter design by picking out the first-Grder transtorna. 2.) 
matrix using Equations 1.6 and 1.11, and then using the 
convolution of coefficients method iteratively on the rows of 
each successive transformation matrix with Equations 1.20 and 
1.21 until the appropriate Ghaurix AG has been generated. 
The next chapter will develop the transformation matrices for 
each of the following transformations: 

- bilinear transformation 

- digital lowpass prototype to digital lowpass filtem 

- @iaqital lowpass prototype to digital highpvass filter 

- digital lowpass prototype to digital bandpass filter 

- diaital lowpass prototype to digital bandstop fiver 

- anaioo lowpass prototype to analog lowpass filter 

- analoo lowpass prototvpe to analoaq highpass filter 

- analog lowpass prototype to analog bandpass filter 

- analoag lowpass prototype to analog bandstop filter. 
We will see that some of these transformations are similiar 
ane consequently. we will Woe have to  aeumee 


Separate transformation matrices for each case. 


Le 


Ii. DEVELOPMENT OF SPECIFIC TRANSFORMATIONS 


A. INTRODUCTION 

The individual transformations are discussed in the 
following sections. The specific form of the transformations 
used in each case are from Reference 2. Since the final 
result will be to develop a computer program to assist in the 
Pett ation Otmes@igital filter coefficients, confining the 
specific transformations to the form in Reference 2, will 
allow the computer program to be used in conjunction with 


Reference 2 aS an ald in designing digital filters. 


B. BILINEAR TRANSFORMATION 

The bilinear transformation may be used to transform an 
eiearod transfer Muerto. Hts), BPemuci = Gigdi1tals) transfer 
femetiwon H(z). Mim clLaistemmataem from H(s) to H(z) is 


described by 


H(z) = H(z)] 
| z-1 (2igAi) 
i 
ei 
Pee see tO BoWetion 1.1, p(x)=z-1 and ag(x)=zt+1l. Since both 
Sf ene transformation polynomials p(x) and q(x) are first 


eraere the badinear transformation wili not change the order 
of the transfer function. Therefore, for an analog transfer 


Mime =Nom of Order mn’, the transformation matrix A(n) will be 


ioe 


(n+1)> vex Betas From Equataon 1.6, the numerator: jen 


denominator polynomial of the transfer™funetron H(z) 125 joaven 


by 

i=n 
a’ (zy= ae’ as *2t+ 7. .+ an '2® =" 2:2 -l] ieee 

ae (255 
where the unprimed coefficients, ai, come from the analog 


transfer function's numerator or denominator polynomial yee 
matrix form, with n=1 and referring “to” Equation “272i 


first order matrix A(i) for the bilinear transformation is 


= cam 
ites 1 | 

A(i) = : (2.28 
-} 1 
= _J 


Usina the convolution of coefficients method developed in 
Chapter 1, the transformation matrix A(nt+1) can be developed 
from the Mac yd Meee ee The results for the bilinear 


transformation matrices A(1l) through A(5) are given in Table 


co 
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TAS iE oy 
Beer TH-ORDERe BILINEAR TRANSFORMATION MATRICES 


Fark S TLaORD EF Ree TO 


aa 


A(1) 


A(2) 


i 
i 


Paes) 





cd cd 


a 
' 


4 


us 
| sm | 


10 


=e 


ca 


oe ee ee | a 


Referring to Table 1 and keeping in mind the Veenvoige 
of coefficients method that produced the matrices, the 
followina relations are apparent and will be useful in 


computer implementation. 


A(n)oi (the first row) = (i) for j=0. 1.2. en (2 ae 
where (4) indicates the binomial coefficient 


A(n)io (the first column) = (-1)* £4£ox 1-073 (2 2 

A(njin (the Last column) ) =a for 1=0,1. 227 (2.6) 

A(n)cn-491 = (-1)"°°-JA(n)i;4 for 1=0,1...,integer(n/2)4+l (eee 

Afn)i:1 = A(n-1),5 + A(m=1)5 (4h Beor 1-0 7 ee (2.8) 
J=L7 23 Re. Oe 


Notice that the values of’ the matrix elements for (the 
biiinear transformation are whole numbers. This is important 
from @ mumerical accuracy standpoint since whole numbers can 
be expressed precisely in floating point form. This means 
that the bilinear transformation will cause a loss of 
accuracy of the original coefficients only because of theme 
multiplications and additions used to caiculate each new 


digital GCoectivei en. 


Ci. LOWPASS TO LOWPASS DIGITAL TRANSFORMATION 
The lowpass to lowpass digital transformation is used to 
transform a lowpaSS prototype digital filter H(z) into a 


lowpass Gigital filter =a The polynomial substitution 
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used 1s 
2) =a Z ) I 
Zo 
j= (Zo) 

ie AZ 
where the subscript p indicates the prototype digital filter 
transfer function. The desiaqn constant a is determined 
[References 2 and 3] by 
sin(@-/2 - 6c '/2) 
Sin(@-/2 - @"/2) (22 On) 
with @- = the prototype critical frequency (usually pi/2) and 
@-'= the desired critical frequency of the lowpass filter. 
Mean, reterrinc to Equation 1.1 with p(x)= z-a and g{z)= 1- 
az, the (feeess digital transformation matrix will not change 
the order of the transfer function and the transformation 
matrix Ain) will be (nt+1) by (nt+1). Using Equations 1.6 and 


[ere aS before, the first order transformation matrix A(1) 


for the lowpass digital transformation is: 


—— ner, 


ec 

At?) = | ela 
|-o 1 
= a 


Once aqdain using the convolution of coefficients method, the 
higher-order transformation matrixes A(n) are developed from 
the lower-order matrix A(n-1). The results are shown in 
Table 2 for the first-order to fifth-order transformation 


matrices. 


a) 


TABLE 2 
FIRST-ORDER TO FIFTB-OKDER UOGWEeAes DiGw am 
TRANSFORMATION MATRICES 


Fe) — ot -a | 
| | 
t é 
a 
A(2)=! 1 -2a a2 | 
| | 
I-a (lta?) -a| 
| | 
| aie i a] 
= wal 
(Se -3a 302 -as 
| =a Cle Zare ) (-2a-a3 ) o | 
| | 
loz C= aa") ee Ores -a | 
| 
| = Bee =3e i | 
= x 
ee -40 G2 ~4q7 a4 | 
| | 
arg (7 gic Za) (~30 =a") (3G “ioade) —3 | 
| | 
mee (~2a-207 (ee ie ee ee a2 | 
oe (Jess oe = Sqr (it see =Oioe 
! | 
ley =e ba? ae ay 
Ref Si=) 4 -5o 1002 -1003 504 ou 
| | 
-9 (1422 |) (-40-602 ) (602 +a4 ) (-40%-a5) a4! 


a? (~2a-3a2) (1+6a24+3a4) (-3a-6a3-a5) (3a2+2a4) -a? | 


| 
—a° (3a4+20%) (—-3@=60%—o5))” (ite eee: (—Za- Sao 7 
ai (~4a%-a5) (602 +404) (-4a-6a3 ) (1#402)  -a | 
| 
a Sag -1Ga3 Lia =e a 
= a= 


18 


Referring to Table 2, the following relations are 


observed and will be used in the computer implementation. 


A(n)oj (the first row) = ( ‘. (~a)J Or T=0F ly ai VAGe sl 2 ) 
A(n)io (the first column) = (-a)! for lon. ...7n (22 13) 
A(n)i 5 =r) < nt) ¢ n= jf) mom, eso, il, . <«.n (2544) 
A(n)i i Seereinia lier yea) tC HO) AH) i C2 sal oe 
Pomii—1, 2, n-l 
Mist 2p es 


Note that the elements OL the lowpass digital 
transformation matrices are not whole numbers as in the 
bilinear transformation case, but are polynomials of the 
design constant a. Numerical error is introduced in this 
Same DY both the ntl additions or multiplications, and the 
error in calculating the matrix elements that multiply the 


Smmainal coefficients. 


I. Meweeoo tO HIGHPASS DIGITAL TRANSFORMATION 
The lowpass to highpass digital transformation is used to 
transform a lowpass prototype to a highpass filter. The 


Bewmymomal Substitution for this.transformation is: 


Hitz) = Ho. (z) Z~-a 
\z=- (216 } 
eZ. 
Since this transformation is the same as that for the lowpass 


to lowpass digital transformation except for the minus sion, 
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we can make the substitution of ‘“-Z) f6r “2 in he origina 
prototype transfer © function and then use the same 
transformation matrices as the lowpass to lowpass case. For 
the lowpass to hiahpass transformation, the design constant a 


is (References 2 and 3) 


cos(@./2 - 6c '/2) 
a= Ss eee (2.099 
cos(@- /2 + @-'/2)} 


with @c - the prototype critical frequency (usually pi/Z)ee 


@c' - the desired critical frequency of the highpass filter. 


Ee LOWPASS TO BANDSTOP DIGITAL TRANSFORMATION 
The polynomial substitution which transforms a Gilqaaae 


lowpass prototype to a bandstop digital filter is (References 


-_ 


Zrangao 
Tio i= 
H(z) = Hertz) 220 = See lS 
i+k i+k (2.18) 
(z= 
20. i-k 
Ll a 2 + eee 
1+k i+k 
cos(O,n '/2 + 8; '/2) 
. = eae ee (2.19) 
cos(@u'/2 - 6); '/2) 
ana 
K = tan(@-. /2)tan(@u '/2-6) '/2) (2.208 
with once aqain, © = the prototype critical frequency 
(usually pi/2)., O.' = the desired upper critical digital 
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Prequency, ae =the desired low critical didqital 


frequency. 


By letting b= 2a/(k+1) ana c=(1-k)/(1i+k). Equatwon 2.18 can 


be simplified to 


H(z) = Hp (2) 


| 
| Zoo DZ + C 
a ee 
t= b2 + ez 
Paain, using Equations 1.6 ana Ald the first-order 


mm@ansformmation matrix 1s 


| 
| 


A(1) = (1 Sales c | 
| | (22) 
| = 1 | 
i i 
Peter a= S2rce the substitution polvnomials are second 


Bane eaten csormatrnior will double the order of the 


ae UTC} geers, tre trensformeaticn 
Meets A’e) *ilt be (n+1) by (2n41). The higher order 
transformation matricies are again developed recursively and 


the results fOr the first-order Xo} third-order 


transformations are given in Table 3. 


a1 


TABLE 3 
FIRST-ORDER™TO THERD-Cl Pee ees fer 
DIGITAL TRANSFORMATION MATRICES 











A(1l) = j}1 =5 e 
Lc = 5 1 
a a 

A(2)= 2 -2b (b2 +2c) -2bc c?2 | 
| | 
e {(-b-bc) (1+b?+c*} 4eb-bem “ec 
lc? -2be (b2 +2c) -2c ae 
L _ 

A(3) = _ -3c (3b2 +3c) (-b3 -6bc) _ 
Ic = «(-2be-b) = s« (1+2b? +be +2c, (=b9=-2b-2bc-2162)) ae 
| 
he (]=2bc-be*)" (eo r2beerZere-) 
l3 -3c2b (3b2 c+3c? ) ee a 
_ : 

(Redundant entries not shown for A(3) due to space 
. ate t sence Entries not shown can be found usSing Equation 
Zo) 


Notice the increased amount of calculation needed to 


determine the matrix elements for the bandstop case. Each 
matrix element is aefunction of 'b' and ‘c’' which are 
themselves functions of the desiaqn constants a and k. Tabie 


3 only gives the transformation matrices for first- to third- 
order. However. to illustrate the amount of calcuJajeaag 


reauired. the matrix element A(5)e25 is 
A(5)en = (-12bc? -6bc? -8b? c-6bc-6b3 -3b-3bc4 -6b?3 c2 -b5 ) (2.2238 


Since each coefficient in the final filter is found by 
Tultiplvino n+l matrix elements by the n+1 origaumes 


Ze 


coefficients and then adding the nt+1 products, the accuracy 
of the design constants and the original coefficients can be 
Critical in the calculation of the final filter coefficients. 


The symmetry relations for the bandstop transformation 


a; 
A(n)io (@ne@frrst Solmmn) = c! HOG 1=O. ...n (Zee a4) 
emg d = Al(N)icn-ir(2n-i1) her Ot oo aN (22 Sy 
a erase 20) 
mera ss = A(n-Mi; - PAW M—-1)ici-1) +°cAenm~-1)i ¢ 1-09 (2.26) 
hore l=, 1,6 .%,n- 1 
J=OR, oe. 2n 


F. LOWPASS TO BANDPASS DIGITAL TRANSFORMATION 
Again, the polynomial substitution which transforms the 


lowpass prototype to a bandpass filter 15 








H(z) = Hp (2) 
2ak Ke 
LE, 
| ean k+1 C2 27a 
iz = - 
Zeek. ke 1] 
ae Zor Ze 
k+1 kee: 
where in this case k = tan(O6/2)cot(@r'/2 - 6:'/2) anda is 
found from the same expression as the bandstop case. With 


eiewcuuolatiedtach of b=Zak/(k+1) a@nd c=(k-1)/(k+1), Equation 


2.2/ becomes 
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H(z) = Hp (z) | 


z* — “Bige + ue 
i (2.20 
1 =—DZ ah ce- 
Since this is the same form as the bandstop 
transformation, we can use the matrices previously developed 


if '-z' is first substituted for 'z' as was done in the 


lowpass and highpass situation. 


G. LOWPASS TO LOWPASS ANALOG TRANSFORMATION 


Tne lowpass to lowpass analog transformation is used to 
transform an analog lowpass prototype to a lowpass analog 
filter. The polynomial substitution is a simple one, namely 

H(s) = Hp (Ss) 
(Qaa2 o>) 
fs S/W 
where we 1s the prewarped critical frequency (we =tan(@- /2) 


and @ is the desired digital critical frequency 


(G-=2(pi)fe/fs). This transformation matrix is found to be 
(1 0 0 0 | 
| 
0 L/w 0° pee : 
A(n) =(¢ 0 1 ee a Oe 
F | 
| . | 
es 
| 
10 0 O 2. . = ii/weem 
oe fen 
cr A(n) = diaaq( 1, l/w .1/Wwe2? : : : »l/we"™?)? throughout. 
(23303) 
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(eee tmaieeeiioethits Case, each new filter coefficient 
Memenads Only on multiplying the old coefficient of s' by 
(we )-' . Therefore, any loss of accuracy due to the 
transformation is due to this one multiplication and the 


Peeor l2mtroduced in calculation of the matrix element (we)-'. 


Ine LOWPASS TO HIGHPASS ANALOG TRANSFORMATION 
The polynomial substitution for the lowpass to hiahvass 


peametormation is 


H(s) = Hp (s) 


c= We /sS Pose) 


t 


Za)  USlTHOmhOUd tiem Smee wecidee t.3]1, the first-order 


mais tormation matrix 1S 


= 
} 


(GDaeS 2: ) 


» 

Ls 

H 

o 
|. Re 


The method of convolution of coefficients can be used to get 
the nicher-order transformation matrices. However, notice 
iia t fetes on 2.32 1s the same matrix as the lowpass to 
lowpass transformation with the reciprocal taken of each non 
zero element, and then each row reversed. This can be shown 


to be the case in general. so that the lowpass to hiaghpass 


transformation matrix is given by 


eS 


| 
J 


° 0 0 1| (2.33) 
| 
1a 0 Ww OM 
ie | 
Ain) ne | 
i. i 
| | 
10 We ha) 0 0 | 
| 
wer 860) 0 0 | 
=! 
Again. each new filter coefficient depends only on one 


multiplication although the calculation of the matrix element 
itself will reauire calculations involving raising the design 
constant. We. to the (n-1)'® power where iis the power of s 


in the final filter transfer Wunctron- 


I. LOWPASS TO BANDPASS ANALOG TRANSFORMATION 
Tne polynomial substitution for transforming a lowpass 


analog prototype to a bandpass analog filter is 


s* + wo? (2.33 


H(s) = Hp (s) | 
| 
ea 


Bs 


where wo is defined as the geometric mean of the passband; 
Wo = (Wuwi)*, and Bis the width of the passband, B=wu-w. 
(wn and wi are the upper and lower critical frequencies 
respectively). Again using Equations 1.6 and 1.11 and 
picking out the appropriate matrix elements, the first-order 


transformation matrix is given by 


Ze 


Wer 


jie Convolution of coefficients technique is again used to 
generate the higher order transformation matrices with the 
meets for the first-order to third-order matrices listed in 
Table 4. As with the bandpass digital transformation, the 
analoc bandpass transformation will double the order of the 
Pmroinal transfer function since the substitution petenonial 
is of second order. 
TABLE 4 


FIRST-ORDER TO THIRD-ORDER BANDPASS 
ANALOG TRANSFORMATION MATRICES 





0 B 0 | 
Af.) = | | 
| Wo ? 0 ca 
= = 
10 0 B2 0 0 | 
| | 
A(z) = |0 Bwo2 0 B 0| 
| 
1Wo 4 0 2Wo * 0 1 | 
Le ) 
coe? 0 0 BS 0 0 0 | 
| | 
0 0 Bema cir a0 B? 0 0 | 
A(3) = | 
° Bwo4 0 2Bwo 2? O B 0} 
| 
lw. 6 0 2we 4 0 3w.-2 0 1 | 
ie ml 
eee eee o tema t™ = 15 cligetly over half empty ened 


that the nonzero elements are sirpnle functions of the desion 


constants we and B. 


oe LOWPASS TO BANDSTOP ANALOG TRANSFORMATION 
For the analoq bandstop transformation, the polynomi am 


Substitution is the reciprocal of the bandpass case: 


Bs (2.36) 


S2 + Wo? 


Usina the same method as before, the first-order 
transformation matrix is found to be: 

lio 2 o 11 (2. aa 

A(l) = | | 

0 B 0 | 

tas _ 
Note that this is the same matrix as the bandpass case but 
with the rows in reverse order. Using the convolution Ga 
coefficients technique the higher-order matrices are 
developed and the results for the first-order to third-order 
transformation matrices are shown in Table 5. Notice that the 
reversal of the order of the rows compared with the bandpass 
case is true in general. Therefore, once again the matrices 


are over half empty and nonzero elements are simple functions 


of B and Wo. 
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eae ae 
Peet On oUrr LO TAI RD-ORDER BANDSTOP 
ANALOG TRANSFORMATION MATKICES 








hw 2 0 Se 
A(1) = | | 
ie B 0! 
{ = 
Wo 4 0 2Wo ¢ 0 ak 
A(2) = |0 Bwot 0 B 0 | 
| | 
0 0 B? 0 0 | 
bai 9 
a ~~] 
Wo 0 2Wo 4 0 3Wo * ¢ 0 1 
O Bwot 0 2Bwo2 0 B 0 
me) = | | 
P 0 B2 Wo 2 0 B2 0 0 
|0 0 0 Be 0 0 
= 


Mee GlelokaAlL OBSERVATIONS 
Now toac  —the@varlous transformation matrices for digital 
filter desian have been developed, it is possible to compare 
the two methods of digital filter desion mentioned earlier 
and depicted in Figure 1. For this comparison, consider an 
analoc lowpass prototype filter. transfer function with 
numerator of order m, and denominator of order n (m less than 
SEreoual torn). 
Pero eowmemeder the digital to digital direct d@sign 
method. Transforminao the analog prototype to a digital 
prototype usino the bilinear transformation will require m+l 


eo 


Mul t1 Dileaear enc and addart Lemme per davei tal prototype 


COcEPeereme tom the numerator polynomial, and Tet a 
Multia i Cations and additiems per daiga tal prototype 
coefficient for the denominator polynomial. Then, each of 


the digital to digital transformations will require 2(n+1) 
multiplications and additions since a result of the bilinear 
transformation is to make n = om. Therefore, going from an 
analog prototype to a final digital filter will regquiame 
(3n+m+4) multiplications andiadda tiens . If lowpass digital 
prototypes are developed prior to the desian process, this 
can be reduced to 2(nt+1) Since then only the digital to 
Gdicital transformations will be used in designing the ame 
falters 

Now consider the analog design method for the same case 
of an mitt -order numerator and n't® -order denominator lowpass 
pretocvoe transicr Bihumetaere For the lowpass and hidqhpass 
case. each new coefficient will require only 1 multiplication 
andwadda tion. For the bandpass and bandstop transformations 
a simple expression for the number of calculations pen 
coefficient is difficult to find, however referring to Tables 
4 ane 5, the number of multiplications and additionsaae 
coefficient is alwavs less than half the order of the 
polynomial to be transformed, on the average. This gives 
approximately (ntn)/2 for the bandpass to bandstop Gaeumem 
transformations. Finally the analog filter must be converted 


to a dicitai filter using the bilinear transformation Whee 


ow 


will require (n+m+2)} @eweulations per coefficient for the 


highpass and lowpass case ame 2 (meemt1) for the bandpass and 


bandstop case. The total number of calculations in using the 
Pmalog desioam method then is®(n+m+4) forthe high/low pass 
and (2.5(nt+m)+2) for the bandpass/stop. These results are 


Mimmearmized mn Table 6. 


TABLE 6 
NUMBER OF MULTIPLICATIONS AND ADDITIONS REQUIRED 
Poe orrr Cire tO FIND THE DIGITAL FILTER 
COEFFICIENTS FROM A LOWPASS ANALOG PROTOTYPE 
(WITH AN M™E ORDER NUMERATOR AND NT#® ORDER DENOMINATOR) 


Analog Design 


high/low pass bandpass/stop 











anaioc transformations : 2 25) (eneem) | 
| 


PemiuumMear transformation 


Roan | n+m+4 feeee O(N my + 2 | 


| 
n+m+2 | 2(n+tm) +2 | 
| 


Digital Design 


heGnytiow’ pass bandpass/stop 





bilinear transformation moe m+ 2 | n+m+2m | 





diqital transformations | Zt) | - 2(n#1) | 


Total 3(n+1)+m+1 | 3(n+1)+m+1 | 


Ba 


Referring to Paple vr, the advantage of one method over 
the other depends upon what type of filter is being designed 
and the order of the prototype transfer functions numerator 
and denominator. Consider designing a “fig th-onden 
low/hiaghpass Chebyshev or Butterworth filter. In this case. 
the order of the numerator polynomial is zero and the analoa 
method will have a slight advantage of 9 multiplications and 
additions versus 12 for the digital design with precalculated 


prototvpves (and 19 for digital desiaqn without precalculated 


prototypes). However, if a bandpass/stop elliptic filter was 
being designed with n=5, then m=4 and the direct digital 
desian method (with precalculated prototypes) has an 


advantage of 12 multiplications and additions compared to 
approximately 24 for the analoa method. 

A final consideration in favor of analog design is the 
complexity of the elements of the transformation mMatmegeeee 
Since the elements of the digital transformation matrauaee 
require much more calculation than the elements of the analog 
transformation matrices or the bilinear transformation 
(compare Tables 1-5), the problem of roundoff error 1s much 
more serious in the case of the digital to digital design 


method. 


ie. IMPLEMENTATION 

The final step is to implement the algorithms developed 
previously in a computer program to assist in the teaching of 
dicital filter design. Since the program is intended to be 


a2 


> - = -_— —_ 
Cae es. SD 


Mm 


Mei eeomm astro etlonal purposes, the numerica. consid 
Giscussead earlier are not considered as important as 
Semeeretal Simplireity. The direct design method was chosen 


for implementation as the design constants can be found 


easily without reference to “prewarping."” The completed 
Prodram. Called “DFCADD' (Digital Filter Computer Aided 
Deyrect —~Design), 1s contained in Appendix A with its own 
documentation. The prooram follows the procedures set forth 
[oeemapter 10 of keference 2. Prototype coefficients are 


Seeeec Wathin the program for first-order to fifth-order 
Butterworth and Chebyshev (1/2,1,&2dB ripples) fli tems. 
While second-order to sixth-order elliptic lowpass prototype 
@fetricienmts are calculated within the program for any 
Mmemewatnectou,or 17/2,1.,2, and 3 €B passbane rippie and -20,- 
eee oe oo, «6aNndlh6h6UCW7O06CUBlhCUsttOpband attenuation, with the 
PeeonOr Nae DrOGra’y developed bv Professor D.E. Kirk of the 
Navai Postoraduate School and References 4 and 5. As a final 
Seaton, the wser May enter his/her own prototype coefficients 
Memmi NerOQea seCMtnm- Order analog prototype. After calculating 
Mmrmatamealetilter coefficients an option exists for the 
Troorar he create a data file, titled “filter”, with the 


Gicgital filter macnitude and phase as a function of the 


{ 


Ceo teMemeEeetency LronezerO tO pi. This data file can then 
Meus Gewet e@nvy Diotting routine to generate aplot of the 


Meecuency response of the filter. 
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APPENDIX 
DFCADD FORTRAN PROGRAM 


DFCADD FORTRAN A 


3630 3 36 3 3 36 3 3 3 IE DC EE 3 38 IC 3 3 IE DC 3E 38 3 3 3 3 3 3 IE 3 30 3 36 3 3 IE 30 36 3 IE 98 3 3 3 30 38 3E 3 3 3 36 3 3 3 3 IE DE 3 3 3E 9 3 3 
% THIS PROGRAM ASSISTS IN THE CALCULATION OF DIGITAL FILTER x 
*% COEFFICIENTS ACCORDING TO THE PROCEDURE SET FORTH IN CHAPTER x 
% 10 CSUMMARIZED ON PG 695) OF - FIRST PRINCIPLES OF DISCRETE x 
% SYSTEMS AND DIGITAL SIGNAL PROCESSING - BY R.D. STRUM AND X& 
* D.E. KIRK CREFERENCEC2) 3} x 
HEH HEE HIE IE DEDEDE 3 JE IED 3 3E EE IE E33 36 3 336 3 3 36 30 36 3 36 30 3 3 30 36 3 3 36 3 SE 36 3 IE OE 36 9 38 3 9 IE 30 36 IE 3E 9 36 3E 3 3 3 
INTEGER ORDER, TYPE, PASS, ORDER2, END, AFORDF 
REAL KPR, LAM, OMEGAC20),V020),A0020),B0(20),B1(20) 
COMPLEX Z,NEWZ, ZN, 2D, HOFZN, HOFZD,HOFZ 
DOUBLE PRECISION TEMP,SUMR 
DOUBLE PRECISION R,50,A01,A02,B01,B02,B11,B12,H0 
DOUBLE PRECISION A03,A04,B05, B04, B15,B14 
DOUBLE PRECISION SFREQ,AFC,DFC,ALPHA,PI,PID4, DFCD2 
DOUBLE PRECISION ALC,AUC, DLC, DUC, 2K,B,C,DLCD2, DUCD2 
DOUBLE PRECISION APROT,BILINR, DPROT,LPHPTR,BPBSTR, DFLTR 
36363 336 3 EE CE JED 3 30 36 3 3 IE 3 30 38 30 38 38 3 38 38 98 30 36 36 38 3 IE 30 38 IE DE 3 38 EE DE 3 E30 38 3 36 30 IE 3 30 3 3 3 3 IE BE IE IE 3 3 3 
FORMAT OF ARRAYS: 
-APROT CONTAINS THE ANALOG PROTOTYPE 
APROTCNORD, PWR) 
-~BILINR IS THE ANALOG-DIGITAL BILINEAR TRANSFORMATION MATRIX 
BILINRCORDER,ROW,COL ) 
-DPROT IS THE DIGITAL PROTOTYPE FROM APROT AND BILINR 
DPROTCNORD, PWR) 
-LPHPTR IS THE LP&HP DIGITAL-~DIGITAL TRANFORMATION MATRIX 
LPHPTRCORDER, ROW,COL) 
-DFLTR IS THE FINAL DIGITAL FILTER 
DFLTRCNORD, PWR) 
-BPBSTR IS THE BP&BS DIGITAL-DIGITAL TRANSFORMATION MATRIX 
BPBSTRCORDER, ROW, COL ) 
O00 B00 88 8800000008808 BB80 IBGE 
DIMENSION APROT(2,0:10) 
DIMENSION BILINR(10,0:10,0:10) 
DIMENSION DPROTC(2,0:10) 
DIMENSION LPHPTRC10,0:10,0:10) 
DIMENSION BPBSTRC10,0:10,0:20) 
DIMENSION DFLTR(2,0:20) 
3636 36 330 36 3636 36 30 36 36 36 36 36 36 336 36 36 336 36 336 36 3036 36 36 36 36 36 336 36 36 36 36 36 3 336 3 3 6 SE 3 3 IE 28 36 BE 38 36 BE 36 38 3 36 38 3 
OPENCUNIT=8,FILE="FILTER DATA',STATUS="UNKNOWN® ) 
20000000000 R8PL P2008 IRB 
PT=3.1415926535898D+00 © 
SPI=3.14159 
PIDG=P174.0D+00 
OOO COCO COCR G80 GC BB CEG 
CLEAR APROT PRIOR TO LOADING 
DO 1 K=1l,2 
DO 2 L=0,10 
APROTCK,L)=0.0D+00 
CONTINUE 
CONTINUE 
FE 3 IEE IEC 3 IEE IE DE IEE 3 E36 3036 36 3 3 336 30 3 36 36 36 36 3 3 36 36 30 36 36 36 36 36 36 36 3 3 36 3 3 36 3E 30 36 36 3 3 36 36 3E 36 36 3 36 3 
HHH KH IE IE DE IE IEE IC IE IE 3 IIE 3 IE BEDE AE 3K 3 IE HE 3 HE HE DE BE 3 IE DE 3 3 IE BE 3 9 IE DE 3 IE 3 3 3 IC 3 3 DE BE BE 3 IC IE 3 IE IE I HE 
BEGIN INTERACTION 
CALL EXCMS C'CLRSCRI') 
VIRITEC6, X) 
WMRITEC6,*) ! 
MRITEC6, X) 
MRITEC6, ®) 
MIRITEC6, ¥) 


THE DIRECT DESIGN METHOD.' 


TIE IE IE IE IE HEH DEI IE DE EE IE IEE IE IE IE IE IEE DE IE IE DEE IE DE DE DE DE DE DE 3E DE DE DE DE DE 3 IE DE 3 3K DE! 
*XWARNING-THIS PROGRAM IS NOT USER FRIENDLY!!! x! 


WRITEC6,%) *% YOU MUST ENTER VALUES OF THE APPROPRIATE ee 
WRITEC6,%*) '* CREAL OR INTEGER), AND PROPER RANGE 

WRITEC6, %) ei ew were xa WRN HNMMIOET RRS 
VIRITEC6, ®) 

WRITEC6,*) * ** ENTER 1 TO CONTINUE xx ' 

READC5,*) LOOK 

CALL EXCMS C'CLRSCRHN*) 

MRITEC6,%*) "WHAT TYPE OF FILTER DO YOU WANT?! 

WRITEC6,*®) ' €1-G ALLOW ONLY UP TO 5TH ORDER PROTOTYPES)' 
VIRITEC6,*) ' 1-BUTTERWNORTH ' 

WRITEC6,%) ' 2-CHEBYSHEV WITH 172 DB RIPPLE ! 

WRITEC6,%*) ' S-CHEBYSHEV WITH 1 DB RIPPLE ! 

WRITEC6,%) ! G-CHEBYSHEV WITH 2 DB RIPPLE’ 
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"THIS IS A PROGRAM TO CALCULATE DIGITAL FILTERS USING! 


NEWO0010 
NEW00020 
NEWO0030 
NEWO00040 
NEWO0050 
NEW00060 
NEWO0070 
NEWO0080 
NEWO00090 
NEW00100 
NEWO00110 
NEWOO120 
NEW00130 
NEW00140 
NEW00150 
NEWOO0160 
NEWO00170 
NEHOO0180 
NEWO0190 
NEWO0200 
NEWO0210 
NEWO0220 
NEW00230 
NEW00240 
NEWOO2Z50 
NENO00260 
NEWO00270 
NEWOO0280 
NEWO0290 
NEW00300 
NEWO0310 
HEWO0320 
NEHO0330 
NEWO00340 
NEW00350 
NEWO00360 
NEHO0370 
NEH00380 
NEHO0390 
NEHW00400 
HEWO0910 
NEM00420 
NEW00430 
NEWO00440 
NEMO00450 
NEHWO0460 
NEHO0470 
NEHO0480 
NEW00490 
HENOO500 
HENW00510 
HENOO520 
NENHOO0530 
NEMO0540 
HEWOO0550 
HEWO0560 
NEWO0570 
NEWO0580 
NEWO0590 
NEHW00600 
NENW00610 
NEHW00620 
NEWO0630 
NEHWO006490 
NEVI00650 
NEVWI00660 
NEW00670 
NEW00680 
NEW00690 
NEW00700 
NEW00710 
NEW00720 


PIE: 


— mC) 
in tN 


DFCADD FORTRAN A 
WRITE(6,%) ! S-ELLIPTIC (ALLOWS UP TO A 6TH ORDER PROTOTYPE) 'NEW00730 
WRITE(6,%) -PASSBAND RIPPLE 0.5,1.0,2.0,0R 3.0 DB! NEW00740 
WRITE(6,%) ' -STOPBAND ATTENTUATION -20,7-30,-40,-50,-60, * NEW00750 
HRITE(6,%) OR -70 DB NEW00760 
WRITE(6,%) ! 6-OTHER - ALLOWS You TO USE UP TO A LOTH ORDER" NEW00770 
WRITE(6,%) ° PROTOTYPE, HOWEVER YOU MUST ENTER THE ANALOG'NEWO00780 
NWRITEC6,%) PROTOTYPE COEFFICIENTS ! NEW00790 
WRITEC6,%) NEW00800 
WRITE(6,%) ! ¥% INTEGER 1-6 xx! NEW00810 
READ(5,%) TYPE NEW00820 
HRM HMMM MMM HHH KKH KKK KK KKK KKK KKK HK RHR HK KR RRR KK RAK KR KR KRKRRRARR NEW00830 
CALL EXCMS ('CLRSCRN') NEW00840 
WRITE(6,%) NEW00850 
WRITE(6,%) "WHAT IS THE ORDER OF THE PROTOTYPE YOU WISH TO USE?" NEW00860 
WRITE(6,%) NEW00870 
WRITE(6,%) ' - MUST BE AN INTEGER IN THE RANGE 1-5 UNLESS YOU' NEW00880 
WRITE(6,%) ! CHOSE TO ENTER YOUR OWN PROTOTYPE (THEN 1-10)"  NEWNO00890 
WRITE(6,%) ° OR YOU CHOSE AN ELLIPTIC FILTER (THEN 1-6)! NEW00900 
WRITE(6,%) * - OR ENTER "0" IF YOU WANT HELP SELECTING THE ORDER' NEW00910 
WRITE(6,%) ° BY CHANGING SPECIFICATION FREQUENCIES TO THE! NEW00920 
WRITE(6,%) ! CORRESPONDING DIGITAL PROTOTYPE FREQUENCIES ° NEW00930 
HRITE(6,%) ! NOTE: YOU WILL NEED PROTOTYPE FREQUENCY '° NEW00940 
HWRITE(6,%) ! RESPONSE PLOTS LIKE THOSE GIVEN IN REF(2) FOR * NEW00950 
HRITE(6,%) ! BUTTERHORTH FILTERS (FIG 10.33, PGS 6918692) ' NEH00960 
WIRITE(6,%) ! AND CHEBYSHEV FILTERS (FIG 10.3%, PGS 696-698) * NEHO00970 
READ(5,%) ORDER NEW00980 
IFCCTYPE.GT.4). AND. (TYPE.LT.8)) THEN NEW00990 
IFCORDER.EQ.1)THEN NE}I01000 
WRITE(6,%) ‘ELLIPTIC FILTERS GIVEN FOR 2ND-9TH ORDER ONLY' NEWO1010 
WRITE(6,%) NEH01020 
GO TO 10 NEW01030 
ENDIF NEH01040 
ENDIF NEH01050 
NEH01060 
CALL EXCMS ('CLRSCRN') NEWO1070 
MRITEC6,%) "DO YOU WANT TO ENTER SPECIFICATION FREQUENCIES* NEH01080 
WRITE(6,%) "AS ANALOG OR DIGITAL FREQUENCIES?! NEN01090 
MRITE(6,%) ! 1-ANALOG ! NEFI01100 
MRITE(6,%) ° 2-DIGITAL ' NEWO1110 
WRITEC6,%) * x INTEGER 1 OR 2 xx! NEH01120 
READ(5,*) AFORDF NEHO1130 
CALL EXCMS ('CLRSCRN*) NEI01140 
IFCAFORDF.EQ.1)THEN NEMO1150 
MRITEC6,%) "WHAT IS THE SAMPLING FREQUENCY IN KHZ?" NEbIO1160 
WRITEC6,*%) *  *% DECIMAL xx! NEM01170 
READ(5,*) SSFREQ NEH01180 
SFREQ=DBLE( SSFREQ) NEH01190 
ENDIF 3 NEH01200 
CALL EXCMS ('CLRSCRN') + NEWO1210 
HRITEC6, ¥) HEHO1220 
WRITEC6,*%) "WHAT TYPE OF PASSBAND DUE YOU WANT?® NEH01230 
MRITEC6,%) ° 1-LOWPASS' NEWO01240 
WRITE(C6,%) 2-HIGHPASS! NEW01250 
WRITEC6,%) *  .3-BANDPASS! NEI01260 
MRITEC6,*) G-BANDSTOP! NEWO1270 
HWRITEC6, X) NEHO1280 
WRITE(6,%) ' xx INTEGER 1-4 xx ! NE-IO01290 
READ(5,%) PASS NEH01300 
FORMULAS FOR HIGHPASS AND LOWPASS NEHO1310 
IFCCPASS.£Q.1).0R. (PASS. EQ.2)) THEN NEW01320 
CALL EXCMS ('CLRSCRN") NEH01330 
IFCAFORDF.EQ.1) THE! NEVIOI 340 
MRITEC6,%) "WHAT IS THE CRITICAL FREQUENCY IN KHZ? ' NEH01350 
MRITE(6,%) ! ~ NEG 3DB PT FOR BUTTERWORTH ' NEI01 360 
WRITEC6,%) ! ~ RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV' NEI01370 
WRITEC6,%) ! xx DECIMAL xx ! NEWO1380 
READ(5,%) SAFC NENO1390 
AFC=DBLE( SAFC) NEHO1400 
DFC=(AFC/SFREQ) XPIX(2.0D+00) NEHO1410 
ENDIF NEVI01420 
IFCAFORDF.EQ.2)THEN NEW01430 
WRITE(6,%) "WHAT IS THE CRITICAL DIGITAL FREQUENCY?' NEW01440 
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Rite: 


ole) 


DFCADD FORTRAN A 


~ NEG 3DB FOR BUTTERWORTH ° 
~ RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV' 
x% DECIMAL XxX ! 


DIRITEC6,%) ° 
WRITEC6,%) ' 
WRITEC6,%) ° 
READ(5,%) SDFC 
DFC=DBLEC(SDFC) 
ENDIF 
DFCD2=DFC/2.0D+00 
CALCULATE THE APPROPRIATE VALUES OF THE DESIGN CONSTANT ALPHA 
FROM TABLE 10.7 REFERENCE (2) AND REFERENCE (3) 
IFCPASS.EQ.1) THEN 
se gue OCDE CDRA ELPS )+CDPCD2) ) 
i 
IFCPASS.EQ.2) THEN 
ee eine DE AHO OGELDS)+CDECD2 ) ) 
NDI 
HE HH HEH HE HE HE ME FE HE HE HE IE HE IE HE SE HE HE HE HE HE HE HE HE HE HE HE HE TE HE HE HK FE I HE HE FE IE EK FE IE SE IE FE FE FE FE HE HE K FE TE HEE EE EE 
OPTION TO SHOM CALCULATED VALUE OF ALPHA 
CALL EXCMS C'CLRSCRN*) 
WRITEC6,%) "DO YOU WISH TO SEE THE CALCULATED VALUE OF ALPHA?’ 
WRITEC6,%) * 8 1-YES/2-NO RR! 
READ(C5,®) LOOK 
IFCLOOK.EQ.1) THEN 
WRITEC6, ¥) 
WRITEC6,903) ALPHA 
WRITEC6,%) *** ENTER 1 TO CONTINUE xx! 
READ(5,%) LOOK 
ENDIF 
ENDIF 
HK HM HH KM MK HM MK HH KM HH MM HE HH HE EE HE HE HE HEE HE KK KK KKK HK KKK KKK KKK 
BANDPASS AND BANDSTOP FORMULAS 
IFCCPASS.EQ.3).0R.CPASS.EQ.4)) THEN 
CALL EXCMS C'CLRSCRH*) 
IFCAFORDF.EQ.1) THEN 
WRITEC6,%) ‘WHAT IS THE LOMIER CRITICAL FREQUENCY IN KHZ?! 
WRITEC6,%) '° - NEG 3DB FOR BUTTERWORTH’ 
WRITEC6,%) ° - RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV ' 
WRITEC6,%) ° ¥x% DECIMAL x! 
READC5,*) SALC 
ALC=DBLECSALC) 
CALL EXCMS C'CLRSCRN!) 
WRITEC6,*%) ‘MHAT IS THE UPPER CRITICAL FREQUENCY IN KHZ?7° 
WRITEC6,%) ' ~ NEG 3DB FOR BUTTERWORTH’ 
WRITEC6,%) ° - RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV ' 
WRITEC6,%) ° R¥ DECIMAL xx! 
READC5,%) SAUC 
AUC=DBLECSAUC) 
DLC=CALC/SFREQ)¥(€2.0D+00)¥PI 
DUC=CAUC/SFREQ)¥(2.0D+00) PI 
ENDIF 
IFCAFORDF.EQ.2) THEN 
WRITEC6,%) ‘WHAT IS THE LOWER CRITICAL DIGITAL FREQUENCY? ' 
WRITEC6,*) ! - NEG 3DB FOR BUTTERWORTH '° 
WRITEC6,%) ° ~ RIPPLE EDGE FOR CHEBYSHEV AND ELLIPTIC’ 
ViRITEC6,%) ° ¥*¥ DECIMAL ®x ° 
READC5,*) SDLC 
DLC=DBLECSDLC) 
CALL EXCMS C'CLRSCRH*) 
WRITEC6,%) * WHAT IS THE UPPER CRITICAL DIGITAL FREQUENCY?! 
WMRITEC6,%) ° ~- HEG 3DB FOR BUTTERWORTH ' 
WRITEC6,%) * - RIPPLE EDGE FOR CHEBYSHEV AND ELLIPTIC’ 
VIRITEC6,*) ° Mx DECIMAL ex ° 
READC5,*) SDUC 
DUC=DBLECSDUC) 
ENDIF 
DLCD2=DLC/2.0D+00 aa 
DUCD2=DUC/2 .0D+00 ‘ 
CALCULATE APPROFRIATE VALUES OF DESIGN CONSTANTS ALPHA AND K 
FROM TABLE 10.7 REFERENCE (2) AND REFERENCE (3) 
ALPHA=DCOS( (DUCD2)+¢(DLCD2) )/DCOS¢ (DUCD2)-(DLCD2) ) 
IFCPASS.EQ.4) THEN 
RK=DTANCPID4) *¥DTANC CDUCD2)-C(DLCD2)) 
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NEWO1450 
NEHO1460 
NEWO1470 
NEHW01480 
NEW014990 
NEWO1500 
NEW01510 
NEHO1520 
NEW01530 
NEW01540 
NEWO1550 
NEWO1560 
NEWO01570 
NEW01580 
NEW01590 
NEW01600 
NEWO1610 
HEWO1620 
NEHO1630 
NEWO1640 
NHEWO1650 
NEWO01660 
HEWO1670 
NEHO1680 
HEMNO01690 
NEWO1700 
HENWO1710 
NEWO1720 
HEWO1730 
NEW01740 
NEbMO1750 
HEMO1760 
NEWO1770 
HEMO1780 
NEWO1790 
NEWO1800 
NEWO1810 
NEWO1820 
HEHO1830 
NEIMO01840 
NEWO1850 
HEWO1860 
NEWO1870 
NEMO1880 
NEHO1890 
HELO1900 
NEMO1910 
HEHWO1920 
HEWO1930 
HEWO1940 
HEMO1950 
HEWO1960 
NEWO1970 
NEWO1980 
HEHO1990 
HEMWO2000 
HEWO2010 
HEMWO2020 
HEWO2030 
HEMO2040 
HEbMO2050 
NEWO2060 
NEMO2070 
NEMO2080 
HEIWO02090 
HEMO2100 
NEr02110 
NEWO2120 
NEWO2130 
NEMO21490 
NEIH02150 
NEHO2160 


FILE: 
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‘eoleleke) 


16 


DFCADD FORTRAN A 
B=((€2.0D+00)XALPHA)/(CRK+1 .0D+00) 
C=((1.0D+00)-RK)/(€1.0D+00)+RK) 

ENDIF 

IFCPASS.EQ.3) THEN 
RK=DTANCPID4)/DTAN( CDUCD2 )-CDLCD2)) 
B=((2 .0D+00 )XALPHAXRK)/CRK+(1.0D+00)) 
NE Oe) ev 


OPTION TO SHOW CALCULATED VALUE OF ALPHA AND K 


CALL EXCMS C'CLRSCRH®) 

WRITEC6,%) "DO YOU WISH TO SEE THE CALCULATED VALUE OF ALPHA ° 
WRITEC6,%) * AND K? ! 

WRITEC6,%) ! RX 1-YES/2-NO RX" 
READ(C5,%*®) LOOK 

IFC(LOOK.EQ.1)THEN 

WRITEC6,%) 

WRITEC6,903) ALPHA 

WRITEC6,90G) RK 

WRITEC6,%) 

WRITEC6,%) '% ENTER 1 TO CONTINUE xx!* 
READ(5,®) LOOK 

ENDIF 


ENDIF 

HH HEHE HEHEHE HE HE HEHE HE HE IEE HEHE DE I HE IE DEH HE HEE HE HE HEI HE HE IEE IE HE IE IE HEE HE HE HE HE DE HE HE HE HE HE HED HE HEH HE HE HEHE HEE EK 
OPTION TO ASSIST IN DETERMINING THE PROTOTYPE ORDER BY CHANGING 
SPECIFICATION FREQUENCIES TO APPROPRIATE PROTOTYPE FREQUENCIES 
FROM EQN 10.274 REFERENCE(]) 

IFCORDER.EQ.O) THEN 


CALL EXCMS C'CLRSCRN®) 
IFCAFORDF.EQ.1) THEN 
WRITEC6,%) "WHAT IS THE ANALOG SPECIFICATION FREQUENCY CKHZ)' 
MRITEC6,%) "THAT YOU WANT TO CONVERT TO A DIGITAL! 
WRITEC6,%) "LOWPASS PROTOTYPE FREQUENCY?! 
WRITEC6,%) * ** DECIMAL Xx! 
READ(5,%) F 
DF=(2.0)*%SPIXCF/SSFREQ) 
ENDIF 
IFCAFORDF.EQ.2)THEN 
WRITEC6,%) "WHAT IS THE DIGITAL SPECIFICATION FREQUENCY ' 
WRITEC6,*) "THAT YOU WAHT TO CONVERT TO A DIGITAL’ 
WRITEC6,*) "LOWPASS PROTOTYPE FREQUENCY?! 
PIRITEC6,*%) *§ %% DECIMAL xX! 
READ(5,%) DF 
ENDIF 
X=COS( DF) 
Y=SINCDF) 
Z2=CMPLXCX,Y) 
IFCPASS.EQ.1) THEN 
AS=REALCALPHA) 
NEWZ=(Z-AS)7001.0)-CAS¥*Z)) 
ENDIF 
IFCPASS.EQ.2) THEN 
AS=REAL CALPHA) 
NEWZ=(-1.0)*%€€Z-AS)7001 .0)-CAS%Z))) 
ENDIF 
IFCPASS.EQ.3) THEN 
BS=REALCB) 
CS=REAL(C) 
Cue ORE K CERI BS REISS) Canaan 
HDI 
IFCPASS.EQ.4) THEN 
BS=REALCB) 
CS=REALC(C) 
me 2 € PSEC FCSN DO ae emcee 
EH 
PTHETA=ABSCATANCAIMAGCNEWZ)/REAL CNEWZ) )) 
IFCPASS.EQ.1)THEN 
IFCF.GT.SAFC) THEN 
PTHETA=SPI-PTHETA 
ENDIF 
ENDIF 
IFCPASS.EQ.2) THEN 


38 


NEW02170 
NEHO2180 
NEW02190 
NEW02200 
NEW02210 
NEW02220 
NEW02230 
NEW02240 
NEW02250 
NEW02260 
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NEW02300 
NEWO2310 
NEW02320 
NEW02330 
NEWO2340 
NEW02350 
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NEW02390 
NEH02400 
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NEWO2460 
NEHW02470 
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NEH02520 
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NEHO2710 
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NEH028 30 
NEH02840 
NEHO2850 
NEW02860 
NEH02870 
NEWO2880 


| ee Se 


eta 


FILE: DFCADD FORTRAN A 


; IFCF.LT. SAFC) THEN 
PTHETA=SPI~PTHETA 
ENDIF 
. ENDIF 
IF(PASS.EQ.3) THEN 
IFCCF.LT.SALC).OR.(F.GT.SAUC)) THEN 
PTHETA=SPI-PTHETA 
ENDIF 
ENDIF 
IFCPASS.EQ.4) THEN 
IFCCF.GT.SALC) AND. (F.LT.SAUC) ) THEN 
PTHETA=SPI-PTHETA 
ENDIF 
ENDIF 
WRITE(6,905) PTHETA 
WRITEC6,%) 
WRITE(6,%) "WANT TO CONVERT ANOTHER FREQUENCY?" 
WRITEC6,%) * XX 1-YES/2-NO Xx! 
READ(5,%) LOOK 
IF(LOOK.EQ.1)GO TO 16 
IFCLOOK.EQ.2) THEN 
WRITEC6,%) 
WRITE(6,%) "WHAT ORDER PROTOTYPE DO YOU WANT TO USE?' 
WRITE(6,%) ' %% INTEGER 1-5 xx ! 
READ(5,%) ORDER 


ENDIF 
ENDIF 
C HH HE HEHE HEHE HEHE HE HEHE HE HE HEHE HE HE HE HEHEHE HE HE HEHEHE HE HEHE HE HEHE HE HEHE HE HE HE HE HE HEHE HE HE HE HE HE HE HE EE HE KE HK KKK 
C LOAD APROT MATRIX WITH APPROPRIATE COEFFICIENTS 
C APROTCTYPE=1) ARE BUTTERWORTH FILTERS 
IFCTYPE.EQ.1) THER 
IFCORDER.EQ.1)THEN 
APROT(1,0)=1.0D+00 
APROT(2,0)=1.0D+00 
APROT(C2,1)=1.0D+00 
ELSEIFCORDER.EQ.2) THEN 
APROT(1,0)=1.0D+00 
APROT(C2,0)=1 .0D+00 
APROTC2,1)=DSQRTC2.0D+00) 
APROT(2,2)=1.0D+00 
ELSEIFCORDER.EQ.3)THEN 
APROT(1,0)=1.0D+00 
APROT(2,0)=1.0D+00 
APROT(C2,1)=2.0D+00 
APROT(2,2)=2.0D+00 
APROT(2,3)=1.0D+00 
ELSEIFCORDER.EQ@.4) THEN 
APROT(1,0)=1.0D+00 
APROTC2,0)=1.0D+00 
APROTC2,1)=2.613126D+00 
APROT(2,2)=3.414214D+00 
APROT(2,3)=2.613126D+00 
APROT(2,4)=1.0D+00 
ELSEIFCORDER.EQ.5) THEN 
APROT(1,0)=1.0D+00 
APROT(2,0)=1.0D+00 
APROT(C2,1)=3.236068D+00 
APROTC2,2)=5.236068D+00 
APROT(2,3)=5.236068D+00 
APROTC2,4)=3.236068D+00 
APROT(2,5)=1.0D+00 
ENDIF 
ENDIF 
C APROTCTYPE=2) ARE CHEBYSHEV FILTERS WITH .5DB RIPPLE 


IFCTYPE.EQ.2) THEN 
IFCORDER.EQ.1) THEN 
APROT(1,0)=2.862775D+00 
APROT(2,0)=2.862775D+00 
APROT(2,1)=1.0D+00 
ELSEIFCORDER.EQ.2) THEN 
APROT(1,0)=1.431388D+00 
APROT(2,0)=1.516203D+00 
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NEWO02890 
NEWO2900 
NEW02910 
NEW02920 
NEW02930 
NEW02940 
NEW02950 
NEW02960 
NEW02970 
NEWO2980 
NEWO02990 
NEWO03000 
NEWO3010 
NEWO3020 
NEWO3030 
NEW03040 
NEW03050 
NEWO3060 
NEW03070 
NEWO03080 
NEW03090 
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NEW03110 
NEWO03120 
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NEWO3140 
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NEWO03160 
NEWO3170 
NEWO03180 
NEW03190 
NE#HO3200 
NEW03210 
NEWO03220 
NEW03230 
NEWO03240 
NEWO3250 
NEW03260 
NEWO3270 
NEbI032380 
NEW03290 
NEHO3300 
NEWO3310 
NEWO3320 
NEWO3330 
HEHWO3340 
NEW03350 
NEW03360 
NEWO3370 
NEH033380 
NEW03390 
NEWO3400 
NEW03410 
NEWO3420 
NEW03430 
HEWO3440 
NEWO3450 
NEbH034960 
NEMO3470 
NEWO34980 
NEW03490 
NEKH03500 
NEWO3510 
NEWO3520 
NEWO3530 
NEWO3540 
NEWO3550 
NEW03560 
NEWO3570 
NEbHO03580 
NEH03590 
NEW03600 


FILE: 


C 


Cc 


DFCADD 


EtlDIF 


APROTCTYPE=3) ARE CHEBYSHEV FILTERS WITH 1DB RIPPLE 


FORTRAN A 


Newt New Ne Ne ee ee ee Pe ee ee te ee FT ee ee te es we TT we we 


to it i) Ot 08 ad et eae a 


Hee OOO 


cc 


D+00 

- 3) THEN 
.715694D+00 
.715694D+00 
-534895D+00 
.252913D+00 
.0D+00 
~EQ.G) THEN 
.357847D+00 
.379051D+00 
-025455D+00 
.716866D+00 
.197386D+00 
.0D+00 
~EQ.5) THEN 
.178923D+00 
.178923D+00 
.752518D+00 
.309575D+00 
.937368D+00 
.172491D+00 
.0D+00 


q 
0 
Q 
7 
1 
5 
2 


——mm OO He 


— eee OO 


IFCTYPE.EQ.3) THEN 
IFCORDER.EQ.1)THEN 
APROT(1,0)=1.196523D+00 
APROTC2,0)=1.196523D+00 
APROT(C2,1)=1.0D+00 
ELSEIFCORDER.EQ.2)THEN 
APROT(1,0)=0.982613D+00 


APROT(C 


APROTCTYPE=49) ARE CHEBYSHEV FILTERS WITH 2DB RIPPLE 


1.102510D+00 
1.097734D+00 
1.0D+00 
~EQ.3) THEN 
0.491307D+00 
0.491307D+00 
1.238409D+00 
0.988341D+00 
1.0D+00 
.EQ.G) THEN 
0.245653D+00 
0.275628D+00 
0.742619D+00 
1.453925D+00 
0.952811D+00 
1.0D+00 
.EQ.5) THEN 
0.122827D+00 
0.122827D+00 
0.580534D+00 
0.974396D+00 
1.688816D+00 
0.936820D+00 
1.0D+00 


IFCTYPE.EQ.4) THEN 
IFCORDER.EQ.1)THEN 
APROT(1,0)=1.307560D+00 
APROT(C2,0)=1.307560D+00 
APROTC2,1)=1.0D+00 
FLSEIFCORDER.EQ.2) THEN 
APROT(C1,0)=0.505803D+00 


APROTC2 
APROTC2 
APROTC2 
ELSEIFCO 
APROT(O] 


, 
, 
, 
R 
, 


02s 


Le 
oa 
DER 
oe 


0.636768D+00 
0.803816D+00 
1.0D+00 
-EQ.3) THEN 
0.326890D+00 


25625D+00. 
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NEW03610 
NEV03620 
NEW03630 
NEW036490 
NEWO03650 
NEN03660 
NEW03670 
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NEW03690 
NEW03700 
NEW03710 
NEW03720 
NEW037 30 
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NEb03790 
NEWO3800 
NEN03810 
NEW03820 
NEWO3830 
NEPNO3840 
NEWO3850 
NEHO3860 
NEWO3870 
NEV03880 
NEM03890 
NEVI03900 
NEHO3910 
NEM03920 
NEHO3930 
NEHW03940 
NENO3950 
NEP03960 
NEH03970 
NEHWO3980 
NEHO3990 
NEHO049000 
NEHWO%9010 
NEH04020 
NEV049030 
NEW04090 
HEW04050 
NEH049060 
NEH04070 
NEWO049080 
NELHOG090 
NEH049100 
NEW04110 
NEb09120 
NEW04130 
NEHO4190 
NEHO4150 
NEHW09160 
NENOG170 
HNEWO04180 
NENWO4190 
NEH04200 
NEH09210 
NEH049220 
NEW049230 
NEW0492490 
NEW049250 
NEW04260 
NEWN09270 
NEbI049280 
NEHW049290 
HEH04¢300 
NEH09310 
NEW04320 


READ(C5,*®) NATTN 
NATTN=-1x*NATTN 
ASUBA=REALCNATTN) 
IFCASUBP.EQ.0.5) THEN 
IFCNATTN.EQ.20) THEN 
IFCORDER.EQ.2)SR=2. 
IFCORDER.EQ.3)SR=1 
IFCORDER.EQ.G9)SR=1 
IFCORDER.EQ.5)5SR=1 
IFCORDER.EQ.6)SR=1 
ENDIF 
IFCNATTIN.EQ.30) THEN 
IFCORDER. EQ. 2)5R=4 
IFCORDER.EQ.3)5R=1 
IFCORDER.EQ.4)SR=1 
IFCORDER.EQ.5)SR=1. 
IFCORDER.EQ.6)SR=1. 
ENDIF 
IFCNATTN.EQ.G0)THEN 
IFCORDER.EQ.2)5R=8 
IFCORDER.EQ.3)SR=2. 
IFCORDER.EQ.4)SR=1 
IFCORDER.EQ.5)5SR=1 
IFCORDER.EQ.6)SR=1 
ENDIF 
IFCNATTN.EQ.50) THEN 


IFCORDER.EQ.2)5R=15. 


IFCORDER.EQ.3)SR=3. 
IFCORDER.EQ.4)5R=2. 
IFCORDER.EQ.5)SR=1 
IFCORDER.EQ.6)SR=1 
ENDIF 
IFCNATTN.EQ.60) THEN 


C FILTERS 
CALCULATED USING A PROGRAM DEVELOPED 


NAVAL POSTGRADUATE SCHOOL, MONTEREY, 


ND RIPPLE DO YOU WANT? CIN DB)’ 
EITHER 0.5, 1.0, 2.0, OR 3.0 Xxx! 


ND ATTENUATION DO YOU WANT? CIN DB)! 
~20,-30,7-40,-50,-60, OR -70 KKK! 
INCLUDE A DECIMAL xxx! 


FILE: DFCADD FORTRAN A 
APROT(2,0)=0.326890D+00 
APROTC2,1)=1.022190D+00 
APROTC2,2)=0.737822D+00 
APROTC2,3)=1.0D+00 

- ELSEIFCORDER.EQ.4)THEN 
APROTC1,0)=0.163945D+00 
APROT(C2,0)=0.205765D+00 
APROTC2,1)=0.516798D+00 
APROTC2,2)=1.2564982D+00 
APROTC2,3)=0.716215D+00 
APROTC2,49)=1.0D+00 

ELSEIFCORDER.EQ.5)THEN 

- APROT(1,0)=0.081723D+00 . 
APROTC2,0)=0.081723D+00 
APROT(2,1)=0.959349D+00 
APROT(2,2)=0.693477D+00 

7 APROTC2,3)=1.4999543D+00 
APROT(2,49)=0.706461D+90 
APROTC2,5)=1 .0D+00 

EWDIF 

- ENDIF 

C APROTCTYPE=5) ARE ELLIPTI 

C ELLIPTIC COEFFICIENTS ARE 

C BY PROFESSOR D. E. KIRK, 

C CALIFORNIA 93943 

A02=0.0D+00 

BO02=0.0D+00 

B12=0.0D+00. 

S0=0.0D+00 

IF ORWPE}E@.5) THEN 
CALL EXCMS C'CLRSCRN!) 
WRITEC6,%) *WHAT PASSBA 
WRITEC6,*) "Xx MUST BE 
READC5,*) ASUBP 
WRITEC6, ®) 
WRITEC6, X) 
WIRITEC6,%) "WHAT STOPBA 
WRITEC6,*®) **** MUST BE 
VIRITEC6,%) *** DO NOT 


76261 


-42189 
-13188 
-09465 
-01555 


.80880 
 Feaoe2 
- 32996 


12912 
05394 


-848925 


71147 


-62842 
.27264 
.12697 


06069 
90430 
06924 


~-48469 
-24101 
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NEWO05030 
NEW05040 
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FILE: 


DFCADD FORTRAN A 


ENDIF 
LFCNATTN.EQ.G0) THEN 
IFCORDER.EQ.2)SR=5.76107 
IFCORDER.EQ.3)SR=2.13923 
IFCORDER.EQ.4)SR=1.40842 
FCORDER.EQ.5)SR=1.16811 
Wee ee 16 
F 
NATTN.EQ@.50) THEN 
FCORDER.EQ.2)5R=10.19181 
FCORDER.EQ.3)SR=3.09137 
FCORDER.EQ@.4)SR=1.75285 
FCORDER.EQ.5)SR=1.33293 
ees EER) 15808 
N 
F 
F 


F 

ATTN.EQ.60) THEN 

CORDER.EQ.2)SR=18 .09398 

CORDER.EQ.3)SR=49.39729 
L1FCORDER.EQ.4)SR=2.299490 
IFCORDER.EQ.5)SR=1.56860 
IFCORDER.EQ.6)SR=1.286953 

ENDIF 

IFCNATTN.EQ.70) THEN 
IFCORDER.EQ.2)SR=32.15951 
LFCORDER.EQ.3)SR=6 .490894 
IFCORDER.EQ.4)SR=2.92363 
1FCORDER.EQ.5)SR=1.88906 
1FCORDER.EQ.6)SR=1 .46330 

ENDIF © 


I 
I 
D 
( 
I 
I 
I 
1 
1 
D 
C 
I 
I 


ENDIF | 
IFCASUBP.EQ.3.0)THEN 


IFCNATT#. EQ, 20) THEN 
IFCORDER.EQ.2)SR=1.73915 
1FCORDER.EQ.3)SR=1.15516 
IFCORDER.EQ.4)SR=1.03853 
IFCORDER.EQ.5)SR=1.00996 
1FCORDER.EQ.6)SR=1.00260 

ENDIF 3 

IFCNATTN.EQ.30) THEN » 
PECORDER.EQ@.2)SR=2.90352 
1FCORDER.EQ.3)SR=1.495814 
IFCORDER.EQ.4)SR=1.19592 
IFCORDER.EQ.5)SR=1.05020 
IFCORDER.EQ.6)5R21.01783 

ENDIF His 

LFCRATTH.EQ.40) THEN 
L1FCORDER.EQ.2)SR=5.05584 
LFCORDER.EQ.3)SR=1.98022 
IFCORDER.EQ.49)SR=1.39663 
IFCORDER.EQ.5)5R=1.13934 
1FCORDER.EQ.6)SR=1.05892 

ENDIF 

LFCNATTN.EQ.50) THEN 
1FCORDER.EQ.2)SR=8.935009 
L1FCORDER.EQ.3)SR=2.79862 
1FCORDER.EQ.4)SR=1.66148 
L1FCORDER.EQ.5)SR=1.28850 
1FCORDER.EQ.6)SR=1.135534 

ENDIF 

IFCHATTN.EQ.60) THEN 
1FCORDER.EQ.2)SR=15.84610 
L1FCORDER.EQ.3)SR=4.03471 
L1FCORDER.EQ.4)SR=2.11596 
LFCORDER.EQ.5)SR=1.50711 
IFCORDER.EQ.6)SR=1.25325 

ENDIF 

L1FCNATTN.EQ.70) THEN . 
1FCORDER.EQ.2)5R=28.15950 
1FCORDER.EQ.3)SR=5.87251 
IFCORDER.EQ.4)SR=2.74749 
1FCORDER.EQ@.5)SR=1.80680 
1FCORDER.EQ.6)SR=1.491799 
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NEHO5840 
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17 


17 


DFCADD 


FORTRAN A 


ENDIF 
ENDIF 
SPI=3.14159 
RK=1.0/SR 
N=ORDER 
EN=N 
KPR=SQRT(1.0-RK¥¥2) 
Q0=0.5*(1.0-SQRTCKPRI)701.04+SQRTCKPR) ) 
Q=Q0+2 .0*Q0%H54+15.0¥Q0%¥9+150 .0XQ0¥X13 
D=(€10.0%*€0.1¥ASUBA)-1.0)7010.0%%(0.1¥ASUBP)-1.0) 
LAM=(1.07(02.0XEN) )¥ALOG((10.0%8(0.05XASUBP)+4+1.0)7€010.0%%(0.05 
XASUBP)-1.0)) 
ISTOP=0 
DENSUM=0 .0 
ENSUM=SINHCLAM) 
M=1 
EM=M 
ENUM=(€(-1 .0) ¥¥M)¥( QHCM¥CM4+1)))¥®SINH(CC2.0XEM+1 .0) LAM) 
DEN=2.0%(C-1.0) ¥¥M) XC QC MHZ) )RXCOSH( 2. OXEMXLAM) 
ENSUM=ENSUM+ENUM 
DENSUM=DENSUM+DEN 
IF(CM.GE.2) THEN 
RN=ABSCENUM/ ENSUM) 
RD=ABS(DEN/ (C1 .0+DENSUM) ) 
IFCCRN.LE.1.0E-8) .AND.CRD.LE.1.0E-8)) THEN 
ISTOP=1 
ENDIF 
ENDIF 
M=M+] 
IFCISTOP.EQ.0)THEN 
GO TO 1? 
ENDIF 
SIGMAO=ABS( (02.0% (Q*0.25)XENSUM)/ (C1 .0+DENSUM) ) 
W=SQRTCC1.0+RKXCSIGMAOX¥2Z) C1 .04+°0S IGMAOX¥2/RK))) 
IN=MODCN, 2) 
IFCIN.EQ.O0)JTHEN 
IR=N/2 
ELSE 
IR=(N-1)72 
ENDIF 
DO 18 I=1,IR 
IFCIN.EQ.O)JTHEN 
EMU=I-0.5 
ELSE 
EMU=I 
ENDIF 
ISTOP=0 
DENHSUM=0.0 
ee ene 
= 
EM=M 
ENUM=((-1.0)%M) ¥( Q¥¥(MHCM4+1]1)))¥SINCC2Z.0XEM+1.0)¥PIXEMU/EN) 
DEN=2.0%(€(-1.0) ¥¥M) ¥( QC M¥¥2))*COS(2.OXEMXPIXEMU/EN) 
ENSUM=ENSUM+ENUM 
DENSUM=DENSUM+DEN 
IFCM.GE.2) THEN 
RN=ABSCENUM/ ENSUE ) 
RD=ABS(DEN/(1.0+DENSUM) ) 
PECCRN.LE.1.0E-8).AND.CRD.LE.1.0E-8)) THEN 
ISTOP=1 
ENDIF 
ENDIF 
M=M+]1] 
IFCISTOP.EQ.0)THEN 
GO TO 19 
ENDIF 
OMEGACI)=(2.0%(Q%*0 .25) XENSUM)/ (1 .0+DENSUM) 
VCID=SQRTCC1 .O-RKXOMEGACI )¥¥2) C1 .0-COMEGACI ) ¥¥2/RK) )) 
AOCI)=1.0/7COMEGACI )¥¥2) 
DN=1.0+CSIGMAOXOMEGACI ) )¥2 
BOCI)=CCSIGMAOXVC I) )¥¥2+COMEGACI ) XW) XX2)7 CDNXZ ) 
B1CIJ=€2.0¥SIGMAOXVCI))/DN 
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Biter 


18 


20 


a 


DFCADD FORTRAN A 
CONTINUE 
HO=1 


0 
IFCIN.EQ.O) THEN 


DO 20 I=1,IR 
HO=HOXBOCI)/AOCI) 

CONTINUE 

HO=HO*X(10.0%*(-0.05*ASUBP) ) 


ELSE 


DO 21 I=1,IR 
HO=HOXBOCI)/A0CI) 

CONTINUE . 

HO=HOXSIGMAO 


ENDIF 


R= 


DBLECSR) 


HO=DBLECHO) 


AO1=DBLECA0(]1 )) 
BOlL=DBLECBO(1)) 
Bl1=DBLECB1(1)) 
IFCCORDER.EQ.2).0R.CORDER.EQ.3)) THEN 
APROT(1,0)=HOXAO] 
APROT(1,1)=0.0D+00 
APROT(1,2)=HO 
APROT(1,3)=0.0D+00 
IFCORDER.EQ.2) THEN 
APROT(2,0)=B01 
APROT(2,1)=B11 
APROTC2,2)=1.0D+00 
ENDIF 


IFCORDER.EQ.3)THEN 
S0=DBLE( SIGMAQ) 
APROT(2,0)=S0*BO1 
APROT(2,1)=B01+(B11%S0) 
APROT(2,2)=S0+B11 
APROT(2,3)=1.0D+00 

ENDIF 


ENDIF 
IFCCORDER.EQ.4).0R.CORDER.EQ.5)) THEN 


AO2=DBLECAQC2)) 

BO2=DBLECBO(2)) 

B12=DBLEC(B1(2)) 

APROT(C1,0)=HOXA01*A02 

APROT(1,1)=0.0D+00 

APROT(1,2)=HO¥CA01+A02) 

APROT(1,3)=0.0D+00 

APROT(1,4)=HO 

APROT(1,5)=0.0D+00 

IFCORDER.EQ.4) THEN 
APROT(2,0)=B01*BO2 
APROTC2,1)=(B11*B02)+(BO01*B12) 
APROT(2,2)=B01+B02+(B11*B12) 
APROTC2,3)=B11+B12 
APROT(2,4)=1.0D+00 

ENDIF 

IFCORDER.EQ.5) THEN 
SO=DBLECSIGMAO ) 
APROT(2,0)=S0*BO1*BO02 
APROTC2,1)=(S0X(CB11*BO2)+(BO01*B12)))+(BO1*B02) 
APROTC2,2)=(S0*(BO01+BO02+(B11*B12) ))+(B11¥BO2)+( BO1¥B12) 
APROTC2,3)=(S0*(B11+B12) )+B01+B02+(B11*B12) 
APROTC2,49)=S0+B11+B12 
APROT(2,5)=1.0D+00 

ENDIF 


ENDIF 


IFCCORDER.EQ.6 


).OR.CORDER.EQ.7))THEN 
AQ2=DBLECAQ(2) 
AOQ3=DBLECA0C3) Wa 
BO2=DBLEC BOC 2) 
BO3=DBLEC BOC 3) 
B12=DBLECB1(2) 
B13=DBLECB1(3) 
APROT(C1,0)=HOXA01 XA02ZKA03 
APROTC1,1)=0.0D+00 


) 
) 
) 
) 
) 
) 


+A * 
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ANMAAAND 


ON ON 
o _ 


DFCADD FORTRAN A 


APROT(1,2)=HOXCCAGZXAD3)+CA01*¥CA02+A03))) 

APROT(1,3)=0.0D+00 

APROTC1,4)=HO*¥(A01+A02+A03) 

APROT(1,5)=0.0D+00 

APROT(1,6)=HO 

IFCCORDER.EQ.6).0R.CORDER.EQ.7)) THEN ° 
APROTC2,0)=B01*BO2*BOS 
APROTC2,1)=CB11¥BO2*B03)+CBO1¥( (BO2*B13)+CB12*BO03))) 
APROTC2,2)=CBO2%B03)+C B11 CC BO2%B13)+CB12%B03)))+ 

(BO1¥((BO2+B03)+(B12*%B13))) 
APROTC2, 3) =BO01*¥(B12+B13)+CB11%C(B02+B03)+(B12*%B13)))+ 
(CBO02*B13)+(B12XB03)) 

APROTC2,4)=B01+(B11¥(B12+B13))+¢€(B02+B03)+¢€B12%B13)) 
APROTC2,5)=B11+B12+B13 
APROT(2,6)=1.0D+00 

ENDIF 

IFCORDER.EQ.7)THEN 
SO=DBLECSIGMAO) 
APROT(2,7)=1.0D+00 
APROTC2,6)=CSOXAPROT (2,6) )+APROTC2, 5) 
APROT(C2,5)=CSOXAPROTC2,5))+APROTC2,49) 
APROTC2,49)=CSOXAPROTC2Z,49) )+APROTC2, 3) 
APROTC2,3)=CSOXAPROTC2,3))4+APROTC2, 2) 
APROTC2,2)=CSOXAPROTC2,2) )+APROTC(2,1) 
APROTC2,1)=CSOXAPROTC2,1))+APROTC(2, 0) 
APROTC2,0)=SOXAPROTC2, 0) 

ENDIF 

ENDIF 
ENDIF 


NORMALIZE THE ELLIPTIC PROTOTYPES 
IFC TYPE.EQ@.5) THEN 
DO 25 I=1,2 
DO 26 J=0,O0ORDER 
APROTCI, J)= =APROTCI, J)/ CODSQRTCR XX J) 
CONTINUE 
CONTINUE 
ENDIF 
APROTCTYPE=6) ALLOWS THE INPUT OF ANY TYPE OF PROTOTYPE FILTER 
BY INPUTING THE COEFFICIENTS OF §S 
IFCTYPE.EQ.6) THEN 
HRITEC6, ¥) 
VIRITEC6,*) "THIS OPTION ALLOWS YOU TO USE ANY ANALOG '! 
WRITEC6,%) * PROTOTYPE OF ORDER 1-10 
WRITEC6,%) *%% COEFFICIENTS MUST BE IN DOUBLE PRECISION: ' 
HRITEC6,%*) ° FORMAT (€0.1234956789D+00) xx ! 
WRITEC6, ¥) 
HRITEC6,*) ‘INPUT THE NUMERATOR COEFFICIENTS FIRST' 
DO 50 I=0,O0ORDER 
WIRITEC6,*®) "WHAT IS THE COEFFICIENT OF Sxx',I,'? 
READC5,*®) APROTC1,I) 
CONTINUE ‘ 
VIRITEC6, ¥) 
CALL EXCMS ('®CLRSCRH') 
PIRITEC6,*) "NOW INPUT THE DENOMINATOR COEFFICIENTS! 
VIRITEC6, *) 
DO 51 J=0,ORDER 
WRITEC6,%) "WHAT IS THE COEFFICIENT OF SxXx',J,'? 
READC5,*®) APROTC2, J) 
CONTINUE 
ENDIF 
HE HE HE HE HEE HE HE HE HE HE DE DE HE HE HE HE HE HE HE HE HE HE HE HE HE DE DE HE DE HE HE HE HE HE HEHE HE HE HE HEE HEHE HE HE HEH HEHE HE EK MMM KX 
HE HE HEH HE HE HE HEH HK HE HE HE HE HEHE HEHE HE HE HE HE HE HEHE HE HE HE HEHE HE HEHE HE HE HE HE HEE HE HE HE HE HE HE HEHE HE HE HEHE HEHE HE HEHEHE HH HK HMR 
DEFINE ELEMENTS OF THE BILINEAR TRANSFORMATION MATRIX 
HE HE HE HE HEH HE HE HE HE HE HE HE HE HE HE HE HE HE HE HEHE HE HE HEHE HEE HE HE HE HE HE HE HE HEH HE HK HH HK HK HH HHH RH HH KKK KKK KKK 
FOR ALL THE BILINEAR TRANSFORATION MATRICES, THE FIRST COLUMN 
HAS ALTERNATING +1'S AND ~1'S. THE LAST COLUMN HAS ALL +#1'°S 
DO 60 I=1,ORDER 
DO 61 J=0,I 
BILINRCI,J,1I)=1.0 
CONTINUE 
CONTINUE 
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FIEE: DFCADD 


OAD 


Q~N 
nn 


OON™N 
tn Oo 


102 


101 
100 


106 
105 


111 


Pi2 


OOO OQ 


FORTRAN A 


FIRST ORDER MATRIX HAS ALREADY BEEN LOADED CONTINUE INTERATION 
UNTIL YOU HAVE THE APPROPRIATE MATRIX 
DO 75 K=2,0RDER 


K oa 
LOAD ELEMENTS IN INTERIOR COLUMNS, ALL ROWS EXCEPT LAST ONE 
DO 76 L=0,KM1 
DO 77 M=1,KM1 
MM1=M-1 
BILINRCK,L,M)=BILINRCKM1,L,M)+BILINRCKM1,L,MM1) 
CONTINUE 
CONTINUE 
NOW LOAD THE LAST ROW 
DO 78 N=1,KMl 
BILINRCK,K,NI=BILINR(K,0,N)¥(0-1.0D+00)%®CK+N)) 
CONTINUE 
CONTINUE 
HE HK HK HH HH HH HK HH HH HH HE HE HK HE HE HE HH HH HH HE HE HE EH HE KE HE HE HE HE HE HE HE HE HE HE HE HE EE EH HR 
THE FOLLOWING CLEARS DPROT PRIOR TO LOADING 
DO 80 K=1,2 
DO 81 L=0,10 
DPROTCK,L)=0.0 
CONTINUE 
CONTINUE 
HE HK HE HH HK HE HE HK HE HH HEH HE KH KK KK HE HE HE HEH HH HE HE HE HK HE HE HEE HE HE KK EE KEE EE EE KKK 
DIGITIZE THE ANALOG PROTOTYPE USING THE BILINEAR TRANSFORMATION 
HEH HE HEH HE HE HE HE HK HE HE HE IE HK HE HE HE KE HE HE HE HE HE HE HE HE EE HE HE HE KE HE HH HE HK HE HEE HH HM KH HH MH HRN 
DO 100 N=1;2 
DO 101 J=0,0RDER 
SUMR=0.0 
DO 102 I=0,0RDER 
TEMP=APROTCN, IO¥BILINRCORDER,I, J) 
SUMR=SUMR+T EMP 
CONTINUE 
DPROTCN, JJ=SUMR 
CONTINUE 
CONTINUE 
HHH HH HH HH HH HH HH HE HE HE KH KH KK HH KK KH HH KK HH HH MH HHH RH KKH 
OPTION TO SHOW DIGITAL PROTOTYPE COEFFICIENTS 
CALL EXCMS C'CLRSCRHN') 
WRITEC6,*%) "THE DIGITAL PROTOTYPE HAS BEEN COMPUTED. DO YOU! 
WRITEC6,%) "WISH TO CHECK THE PROTOTYPE COEFFICIENTS?! 
Dies XX L-YES/2-NO XX‘ 
READC5,*®) LOOK 
IFCLOOK.E£Q.1) THEN 
MAKE DENOMINATOR COEFFICIENT OF Z®XORDER = 1.0 
DO 105 K=1,2 
DO 106 L=0,0RDER 
DPROTCK,LJ=DPROTCK,LI/DPROTC2,0RDER) 
CONTINUE 
CONTINUE 
CALL EXCMS C'CLRSCRN'!) 
WRITEC6,*®) "THE DIGITAL PROTOTYPE NUMERATOR COEFFICIENTS ARE:' 
DO 111 N=ORDER,0O,-1 
WRITEC6,902) DPROTC1,N),N 
CONTINUE 
WRITEC6, ®) 
WMRITEC6,*®) "THE DIGITAL PROTOTYPE DENOMINATOR COEFFICIENTS ARE:' 
DO 112 M=ORDER,O,-1 
WRITEC6,902) DPROT(2,M),M 
CONTINUE 
WRITEC6, ®) 
WRITEC6,®) "® 1 TO CONTINUE x! 
READ(C5,*®) LOOK 
ENDIF 
HE HHH HE HE HEH HH HE HE HE HE HE HE HE HE DE HE HE HE HEH HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HH EH RH 
NOW DO THE DIGITAL TO DIGITAL TRANSFORMATIONS 
HEH HE HH HH HEH HE HHH HEH HH HE HEH HH HK HEH HE HEH HK HE HH HEHE HH HHH HHH HH HH HH HH HHH HHH HH HH HH HH 
IFCCPASS.£EQ.1).0R.CPASS.EQ.2)) THEN 
DEFINE ELEMENTS OF THE LOWPASS AND HIGHPASS DIGITAL-DIGITAL 
TRANSFORMATION MATRIX 
DEFINE THE ELEMENTS FOR FIRST AND LAST COLUMNS 
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181 
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186 
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O09 aM 


226 
225 


DFCADD FORTRAN A 


DO 170 I=1,10 
DO 171 J=0,1 
LPHPTRCI,J,0)=(€((-1.0)*ALPHA)®XJ) 
IMJ=I-J : 
LPHPTRCI,IMJ,I)=LPHPTRiI,J,0) 
CONTINUE 
CONTINUE 
FIRST ORDER MATRIX HAS 2FEN DEFINED ABOVE 
SECOND TO TENTH ORDER BELOW 
DO 175 K=2,10 
KM1=K-1 
DO 176 L=0,KM1 
DO 177 M=1,KM1 
MM1=M-1 


LPHPTRCK,L,M)=LPHPTRCKM1,L,M)+CCLPHPTRCKM1,L,MM1))*C-1.0D+00) 
C¥ALPHA) 


CONTINUE 
CONTINUE 
DO 178 N=1,KM1 
KMN=K-N 
LPHPTRCK,K,N)=LPHPTRCK,0,KMN) 
CONTINUE 
CONTINUE 
HE HH HE HE HE HE HE HK HE HE HE HE HE HE HE DE HE DE DE DE HE HE HE IE HK HE DE DE HE HE DE HE DE HE HE HE DE UE EE EE HE EEE HE EE EE EEE OE EE RH 
CALCULATE THE FINAL FILTER 
FOR HIGHPASS, THE PROTOTYPE MUST SUBSTITUTE (-Z) FOR (2) 
IFCPASS.EQ.2) THEN 
DO 180 K=1,2 
DO 181 L=0,0RDER 
DPROTCK,L)=DPROTCK,LIXCC-1.0)%XL) 
CONTINUE 
CONTINUE 
ENDIF 
DO 185 K=1,2 
DO 186 L=0,0ORDER 
SUMR=0.0D+00 
DO 187 M=0,ORDER 
TEMP=DPROTCK,M)XLPHPTRCORDER,M,L) 
SUMR=SUMR+TEMP 
CONTINUE 
DFLTRCK,LI=SUMR 
CONTINUE 
CONTINUE 
ENDIF 
HE HK KH HK HE HEHE HEE HE DE HE HE HE DE HE DE DE DE HE DE DE DEE HE HE HE DE DE HE HE EI HEE HEE HE EEE HE HE HEE HEE EE HE EH EK EK 
BANDPASS AND BANDSTOP FORMULAS 
IFCRRASS.EQ.3).0R.CPASS.EQ@ S07 RiEN 
LOAD THE BANDPASS/BANDSTOP TRANSFORMATION MATRIX 
FIRST LOAD THE FIRST AND LAST COLUMNS 
DO 225 I=1,10 
DO 226 J=0,]1 
BPBSTRCI,J, 0) =CxxJ 
I2=2xI 
IMJ=I-J 
BPBSTRCI,IMJ,1I2)=BPBSTR(1I,J,0) 
CONTINUE 
CONTINUE 
LOAD THE REST OF THE FIRST ORDER MATRIX 
BPBSTRC1,0,1)=(-1.0D+00)*B 
BPBSTRC1,1,1)=(-1.0D+00)XB 
LOAD THE 2ND-10TH ORDER TRANSFORMATION MATRICES 
DO 230 K=2,10 
K2=2*K 
K2M1=Ke-1 
KM1=K-1 
K2M2=K2-2 
KM12=2*KM1 
KM12M1=KMi2-1 
LOAD THE SECOND AND SECOND TO LAST COLUMNS 
DO 231 L=0,KM1 
BPBSTRCK,L,1)=BPBSTRCKM1,L,1)+CBPBSTRCKM1,L,0)*(-1.0D+00)*B) 
BPBSTR(K,0,K2M1)=CBPBSTRCKM1,0,KM12)%C-1.0D+00)XB)+ 


48 


NEW09370 
NEWO09380 
NEW09390 
NEW09400 
NEW09410 
NEW09420 
NEW09430 
NEW09440 
NEW09450 
NEW09460 
NEW094970 
NEW09480 
NEW09490 
NEW09500 
NEW09510 
NEWO9520 
NEW09530 
NEW09540 
NEWO09550 
NEWO09560 
NEVIO9570 
NEWO09580 
NEWO9590 
NEW09600 
NEW09610 
NEV09620 
NEV09630 
NEW09640 
NEHWO9650 
NEW09660 
NEW09670 
NEHW09680 
NEW09690 
NEWO9700 
NEWO9710 
NEWO9720 
NEW09730 
NEWO09740 
NEWO9750 
NEHO9760 
NEWO09770 
NEWO9780 
NEWO9790 
NEWO09800 
NEWO09810 
NEWHO9820 
NEHWO9830 
NEVHl098490 
NEWO09850 
NEWMO09860 
NEVWO9870 
NEWO9880 
NEVHlo9890 
NEW09900 
NEVH09910 
NEWO9920 
NEWO9930 
NEW09940 
NEWO9950 
NEH09960 
NEWO9970 
NEW09980 
NEW09990 
NEW10000 
NEW10010 
NEW10020 
NEW10030 
NEW10040 
NEW10050 
NEW10060 
NEW10070 
NEW10080 


FILE: 


Zl 
C 


250 


247 


246 
245 


251 


250 


300 


301 


© 


DFCADD FORTRAN A 


ee Oe ct xe) 
BPBSTRCK,KML,K2M1)=BPBSTR(K,L,1) 
CONTINUE 
NOW LOAD INTERIOR ELEMENTS EXCEPT FOR LAST ROW 
DO 236 L=0,KM1 
DO 237 M=2,K2M2 
MM1=M-1 
MM2=M-2 
BPBSTRCK,L,M)=BPBSTRCKM1,L,M)+(BPBSTRCKM1,L,MM1) ®(-1 .0D+00) 
BJ+(BPBSTRCKM1,L,MM2Z) ®C) 
CONTINUE 
CONTINUE 
NOW FILL THE LAST ROW 
DO 238 N=1,K2M1 


K2MN=K2-N ‘ 

BPBSTRCK,K,K2MN)=BPBSTR(CK,0,N) 
CONTINUE 
CONTINUE 


HEH HE HE HE HE HE HE HE IE HE HE HE HE HE HE HE HE HE IE HE IE HE IE HE HE HE HE IE HE IE HE HE HE HE HE HE HE HE IE HE IE HE HE HE HE HE HEE HEE HEE HE EE HE EE KKK 
CALCULATE THE FINAL FIt £ER?: 
FOR THE BANDPASS CASE REPLACE (2) WITH (-2Z) 
IFCPASS.EQ.3) THEN 
DO 240 M=1,2 
DO 241 N=0,Q0RDER 
DPROT(M,N)=DPROTCM,N) ®((-1 .0D+00 ) X¥N) 
CONTINUE 
CONTINUE: 
ENDIF ° 
ORDERZ=2*ORDER 
DO 245 K=1,2 
DO 246 L=0,0RDERZ2 
SUMR=0.0D+00 
DO 247 M=0,0RDER 
TEMP=DPROT(CK,M)*BPBSTRCORDER,M,L) 
SUMR=SUMR+TEMP 
CONTINUE 
DFLTR(K,L)=SUMR 
CONTINUE 
CONTINUE 
ENDIF 
HM HMH HH HHH HHH HH KK 
OUTPUT OPTIONS 
KKK K MK MH K KH KKH KK 
IFCCPASS.EQ.3).0R.CPASS.EQ.9)) THEN 
ENUD=2*ORDER 
FNUDIF 
IFCCPASS.EQ.1).0R.CPASS.EQ.2)) THEN 
END=ORDER 
ENDIF 
MAKE DENOMINATOR COEFFICIENT OF Z¥®XEND EQUAL TO 1 
DO 250 M=1,2 
DO 251 N=0,END 
DFLTROM,N)J=DFLTROM,ND/DFLTR(2Z, END) 
CONTINUE 
CONTINUE 
CALL EXCMS C'CLRSCRN'*) 


WRITEC6, ®) 
"THE FINAL DIGITAL FILTER IS:'* 


WRITEC6, ¥) 
WRITEC6, *) 

VIRITEC6,%) * NUMERATOR COEFFICIENTS: '! 
DO 300 I=END,0,-1 

WIRITEC6,902) DFLTR(1,1),1 

CONTINUE 
WRITEC6,%) * DENOMINATOR COEFFICIENTS: ° 

DO 301 J=END,0,-1 

WRITEC6,902) DFLTR(2,J),J 

CONTINUE 

MRITEC6,*) 

WRITEC6,*®) * ®% ENTER 1 TO CONTINUE xx! 

READ(5,*®) LOOK 

HEHE IE IE IE FEI IE IE IE IESE IEE 3 JE 3 3 3 3 IEE 3 IE DE 3 IE DE I BC 3 I 3 I IE BE I BC 3 IE BC IE BE 3 BE I BE 3 9 EBC HE BE IE IE BE I BE IE IE BE IE 3 


49 


NEW10090 
NEW10100 
NEW10110 
NEW10120 
NEW10130 


NEW10220 
NEW10230 
NEW10240 
NEW10250 
NEW10260 
NEW10270 
NEW10280 
NEW10290 
NEVW10300 
NEW10310 
NEW10320 
NEW10330 
NEW10340 
NEW10350 
NEW10360 
NEW10370 
NEW10380 
NEW10390 
NEW10400 
NEW10410 
NEW10420 
NEW10430 
NEW1049490 
NEW10450 
NEW10460 
NEW104970 
NEW10480 
NEW104990 
NEW10500 
NEW10510 
NEW10520 
NEW10530 
NEW10540 
NEW10550 
NEW10560 
NEIW10570 
NEW10580 
NEW10590 
NEW10600 
HEW10610 
NEW10620 
NEW10630 
NEW10640 
NEW10650 
NEW10660 
NEW10670 
NEW10680 
NEW10690 
NEW10700 
NEW10710 
NEW10720 
NEW10730 
NEW10740 
NEW10750 
NEW10760 
NEW10770 
NEW10780 
NEW10790 
NEW10800 


Pies 


C 


310 


320 


DFCADD FORTRAN A 


OPTION TO SPOT CHECK FILTER OUTPUT MAGNITUDE AND FREQUENCY NEW10810 
CALL EXCMS C*CLRSCRN'*) NEW10820 
WRITEC6,*%) "THE TRANSFER FUNCTION FOR THIS FILTER HAS BEEN' NEW10830 
WRITEC6,%) "COMPUTED. DO YOU WANT TO SPOT CHECK THE MAGNITUDE' NEW10840 
WRITEC6,%) "AND PHASE OF THE FILTER OUTPUT FOR A GIVEN °* NEW10850 
WRITEC6,%) *FREQUENCY?'° NEW10860 
WRITEC6,%) * XX I-YES/2-NO xx * NEW10870 
READ(C5,*) LOOK NEW10880 
IFCLOOK.EQ.1) THEN NEW10890 
CALL EXCMS CfCLRSCRN') NEW10900 
IFCAFORDF .EQ.1) THEN NEW10910 
WRITEC6,%) "WHAT ANALOG FREQUENCYCIN KHZ)' NEW10920 
WRITEC6,%) * DO YOU WANT TO CHECK?! NEW10930 
WRITEC6,%) * *xX% DECIMAL xx° NEW10940 
READC5,*®) F NEW10950 
DF=(F/SSFREQ) *2.0SPI NEW10960 
ENDIF NEW10970 
IF CAFORDF.EQ.2) THEN NEW10980 
WRITEC6,%) "WHAT DIGITAL FREQUENCY DO YOU WANT TO CHECKT'* NEW10990 
WRITEC6,%) * *X DECIMAL Xx * NEW11000 
READ(C5,%) DF NEW11010 
ENDIF NEW11020 
X=COSC DF) NEW11030 
Y=SINCDF) NEW11040 
Z=CMPLXCX,Y) NEW11050 
XN=REALCDFLTR(1,0)) NEW11060 
YH=0.0 NEW11070 
XD=REALCDFLTR(2,0)) NEW11080 
YD=0.0 NEW11090 
HOFZN=CMPLXCXN, YN) NEW11100 
HOFZD=CMPLXCXD, YD) NEW11110 
DO 320 N=1,END NEW11120 
ZN=REALCDFLTRC1,N) 2 ¥C ZN) NEW11130 
ZD=REALCDFLTR(2,H))¥C ZN) NEW11140 
HOFZN=HOFZN+ZN NEW11150 
HOFZD=HOFZD+ZD NEW11160 
CONTINUE NEW11170 
HOFZ=HOFZN/HOFZD NEHW11180 
FMAG=SQRTCCREAL CHOFZ)*®*2)+CAIMAGCHOFZ)XxX2)) NEW11190 
FMAGDB=20.0QXALOGIOQCFMAG) NEW11200 
FPHSE=ATANCAIMAGCHOFZ)/REAL CHOFZ)) NEW11210 
WRITEC6,%) "THE FILTER OUTPUT WILL BE:' NEbW11220 
WRITEC6,919) FMAGDB NEW11230 
WRITEC6,915) FPHSE NEW112490 
WRITEC6, ®) NEW11250 
WRITEC6,*) *DO YOU WANT TO CHECK ANOTHER FREQUENCY?! NEW11260 
WRITEC6,*) * XX L-YES/2-NO xx ° NEW11270 
READC5,*) LOOK NEMW11280 
IFCLOOK.EQ.1)G0 TO 310 NEW11290 
ENDIF NEW11300 
HE HEH HEH HE HEHE HE HEHE HEHEHE HE HE HE HE HEHEHE DE DE HE HE HE HEHE HE DE DE HE IE HE HE HE DE DE IE DE DE DE HEHE DE DE HE DEE DEE EEE HHH HRKK NEYI1I310 
OPTION TO LOAD FREQUENCY RESPONSE INTO A DATA FILE NEHW11320 
CALL EXCMS CTCLRSCRH') NEW11330 
WRITEC6,*) "DO YOU WANT THE FILTERS FREQUENCY RESPONSE ! NEW11340 
WRITEC6,%) "DATA TO BE ENTERED INTO A DATA FILE FOR USE‘ NEW11350 
WRITEC6,*) "BY A PLOTTING PROGRAM? ! NEW11360 
WRITEC6,*) * NOTE: THE FILE NAME WILL BE FILTER, FILE TYPE! NEW11370 
WRITEC6,%) ° WILL BE DATA, 100 POINTS WILL BE CALCULATED! NEMW11380 
WRITEC6,%) * FOR A DIGITAL FREQUENCY FROM O TO PI, THE FIRST* NEW11390 
WRITEC6,%) ° COLUMN WILL BE THE DIGITAL FREQUENCY CEXPRESSED ‘NEW114900 
WRITEC6,x%) * AS A FRACTION OF THE SAMPLING FREQUENCY 0-172) ' NEW11410 
WRITEC6,%) THE SECOND COLUMN WILL BE THE MAGNITUDE, AND‘ NEW11420 
WRITEC6,%) ° THE THIRD COLUMN WILL BE THE PHASEC DEGREES)! NEW11430 
WRITEC6,%) ° *¥*% 1-MAKE DATA FILE/2-DON''T BOTHER xXx! - NEW11440 
READ(C5,*®) LOOK NEW11450 
IFCLOOK.EQ.1) THEN NEW11460 
DO 350 I=0,100 NEW11470 
DF=REALCI)*¢€5.0/71000.0) NEW11480 
DFR=REALCI)X¢(SPI4100.0) NEW114990 
X=COSCDFR) NEW11500 
Y=SINCDFR) =o NEW11510 
Z2=CMPLX(X,Y) ‘ NEW11520 
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FILE: 


S51 


350 


901 
902 
903 
904 
905 
914 
915 
916 


DFCADD FORTRAN A 


oN eee Re »0)) 
XD=REAL CDFLTR(2,0)) 
YD=0.0 


HOFZN=CMPLXCXN, YN) 
HOFZD=CMPLX(XD,YD) 
DO 351 N=1,END 
ZN=REALCDFLTRC1,N) ) (ZN) 
ZD=REALCDFLTR(2,N) I XC ZRKN) 
HOFZN=HOFZN+ZN 
HOFZD=HOFZD+Z2D 
CONTINUE 
HOFZ=HOFZN/HOFZD 
FMAG=SQRTC CREAL CHOFZ )**2 )+CAIMAGCHOFZ) ®*2) ) 
IFCREALCHOFZ).NE.0.0)THEN 
FPHSER=ATANCAIMAGCHOFZ)/REALCHOFZ)) 
ELSEIFCREALCHOFZ).EQ.0.0)THEN 
IFCAIMAGCHOFZ).GT.0.0)FPHSER=1.570796 
IFCAIMAGCHOFZ).LT.0.0)FPHSER=-1.570796 
IFCAIMAGCHOFZ).EQ.0.0)FPHSER=0.0 
ENDIF 
FPHSE=FPHSER¥G180.0/SPI) 
WRITEC8,916) DF, FMAG, FPHSE 
CONTINUE 
ENDIF 
HE FE HE HEHE HEHE HE IE HE HE FE FE HE HE IE FE HE HE IE HE FE IE FE HE FE HE FE FE HE HE HE HE HE HE HE FE HE HE FE HE HE HE HE HE HE IE FE HE IE HE HE HE IE HE HE HEHE HE HE HHH 
FORMAT(C5X,F10.5,° Sxx',I1) 
FORMAT(5X,F10.5,' Zxx",11) 
FORMAT( 3X, 'ALPHA= ',F9.4) 
FORMATC3X,'K= ',F9.9) 
FORMAT(C3X,'THE LOWPASS DIGITAL PROTOTYPE FREQ. IS 
FORMATC3X,'"MAGNITUDE = ',F8&.2,' DB') 
FORHATC3X,*PHASE = ',F6.2," RADIANS") 
FORMAT(3X,F4.3,10X,F5.3,10X,F6.2) 
STOP 
END 


Ve Oo) 


ol 


NEW11530 
NEW115490 
NEW11550 
NEW11560 
NEW11570 
NEW11580 
NEW11590 
NEW11600 
NEW11610 
NEW11620 
NEW11630 
NEW11640 
NEW11650 
NEW11660 
NEW11670 
NEW11680 
NEW11690 
NEW11700 
NEW11710 
NEW11720 
NEW11730 
NEW11740 
NEW11750 
NEW11760 
NEW11770 
NEW11780 
NEW11790 
NEW11800 
NEW11810 
NEW11820 
NEW11830 
NEW11890 
NEW11850 
NEW11860 
NEW11870 
NEW11880 
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